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SECTION  I 


INTRODUCTION 


The  guidance  and  control  field  has  traditionally  focused  on  continu¬ 
ous  or  analog  control  systems  represented  in  the  Laplace  or  s-domain  or 
in  a  state-space  model.  Today,  the  increasing  use  and  popularity  of  the 
digital  computer  as  a  system  component  has  provided  the  major  impetus  to 
the  theoretical  as  veil  as  the  practical  interest  in  sampled-data  or 
discrete  control  systems.  The  basic  problem  facing  the  control  engineer 
is  obtaining  valid  discrete  models  of  complex,  closed-loop  hybrid 
systems  (L.e.,  systems  containing  both  analog  and  discrete  elements). 
These  discrete  models  must  be  in  a  convenient  form  that  can  be  readily 
analyzed  using  the  analysis  and  synthesis  tools  available  today.  Valid 
discrete  models  for  hybrid  systems  can  be  obtained  using  the  z- ,  w-,  or 
w'-transform.  The  z-transform  is  a  logical  extension  of  the  Laplace 
transform  and  can  be  used  to  handle  sampled-data  systems.  The  w-  and 
w'-transforms  are  related  to  the  z-transform  through  simple  bilinear 
transformations.  Discrete  models  expressed  in  the  z-,  w-,  or  w'-plane 
define  the  continuous  variables  in  a  hybrid  system  at  the  sampling 
instants  and  completely  describe  the  inherently  discrete  variables  asso¬ 
ciated  with  digital  elements. 

The  two  computer  programs  presented  in  this  volume  provide  some  of 
the  basic  digital  implementation  tools  required  in  the  analysis  and  syn¬ 
thesis  of  hybrid  systems.  The  DISCRET  computer  program  converts  a 
general  analog  or  continuous  model  expressed  in  the  s-plar.e  to  the  z-, 
w-,  or  w'-plane.  DISCRET  can  calculate  the  standard,  de  ayed,  or  ad¬ 
vanced  discrete  transform.  Data  holds  including  the  zero  order,  first 
order,  second  order,  and  slewer  can  be  inserted  into  the  transforma¬ 
tion.  The  second  program,  TXCONV,  implements  the  conversion  of  a  high- 
rate  discrete  transform  to  a  low-rate  discrete  transform.  The  general 
input/output  structure  of  these  two  computer  programs  is  shown  in 
Fig.  1.  In  DISCRET,  the  input  is  an  s-plane  transfer  function  and  the 
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Figure  1.  General  Structure  of  DISCRET  and  TXCONV 
Computer  Programs 
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output  Is  a  selected  discrete  transfer  function.  For  TXCONV,  the  input 
Is  a  high-rate  discrete  transfer  function  and  the  output  a  low-rate  dis¬ 
crete  transfer  function  in  the  z-,  w- ,  w'-plane. 

In  Section  II,  a  review  of  basic  sampled-data  theory  is  presented. 
This  section  provides  the  necessary  background  information  for  succeed¬ 
ing  sections.  The  mathematical  development  is  based  on  the  assumption 
that  the  sampling  process  can  be  described  as  the  amplitude  modulation 
of  an  impulse  train  by  the  input  signal.  This  assumption  greatly  sim¬ 
plifies  sampled-data  theory  and  is  valid  for  most  practical  engineering 
applications. 

The  transform  conversion  expressions  mechanized  in  the  TXCONV  compu¬ 
ter  program  are  developed  in  Section  III.  These  expressions  allow  a 
high-rate  transform  in  the  z-,  w-,  or  w'-plane  to  be  converted  to  a  low- 
rate  discrete  transform.  The  fundamental  definition  of  the  z-transform 
and  the  discrete  inversion  integral  developed  in  Section  II  form  the 
basis  for  the  analytical  development.  The  computer  mechanization  of  the 
transform  conversion  is  based  on  the  practical  calculation  of  the 
residues  of  a  complex  integral.  The  residues  for  this  integral  are 
calculated  in  an  unconventional  manner  using  a  limiting  process  via 
L'Hopital's  rule.  This  method  simplifies  the  mechanization  scheme  and 
leads  to  a  closed-form  solution. 

Sections  IV  and  V  provide  a  detailed  description  of  the  DISCRET  and 
TXCONV  computer  programs,  respectively.  This  includes  the  available 
program  options,  theoretical  basis  for  the  mechanization  algorithms, 
general  and  detailed  program  structure,  required  input  data,  and  typical 
program  output.  Source  listings  for  these  two  computer  programs  are 
contained  in  Volume  III. 
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SECTION  II 


REVIEW  OF  FUNDAMENTAL  SAMPLED-DATA  THEORY 

A.  INTRODUCTION 

A  quick  review  of  the  fundamental  principles  of  sampled-data  theory 
is  presented  in  this  section  (Refs.  1—10) .  This  background  information 
will  be  used  in  succeeding  sections  to  develop  the  analytical  expres¬ 
sions  mechanized  in  the  DISCRET  and  TXCONV  computer  programs.  The  basic 
theory  for  sampled-data  or  discrete  systems  vas  developed  over  20  years 
ago  and  remains  intact  today.  Practical  sampled-data  theory  is  based  on 
the  assumption  that  t.  actual  sampling  operation  can  be  modeled  as  the 
amplitude  modulation  of  an  impulse  train.  This  central  concept  greatly 
simplifies  the  analysis  and  synthesis  of  sampled-data  systems.  For¬ 
tunately,  this  view  of  the  sampling  process  is  valid  for  roost  practical 
systems  and  use  of  this  theory  is  normally  considered  exacts 

Sampled-data  systems  generally  contain  both  continuous  and  discrete 
elements.  The  z-transform  provides  a  unified  analysis  and  synthesis 
technique  for  these  hybrid  systems.  For  a  sampled  continuous  element, 
the  z-transform  can  be  considered  as  the  Laplace  transform  of  an  impulse 
sequence  (impulse  train)  where  the  area  or  strength  of  the  individual 
impulses  equal  the  value  of  the  continuous  time  function  at  each  dis¬ 
crete  sampling  instant.  An  alternate  viewpoint  is  to  consider  the 
exponent  in  the  z~n  delay  operator  as  an  ordering  variable  for  a  number 
sequence  (or  a  sequence  of  discrete  signal  values)  where  the  coefficient 
for  the  z~n  terms  equals  the  value  of  the  number  sequence  (or  the.  dis¬ 
crete  signal)  at  the  nth  discrete  time  instance.  This  viewpoint  allows 
the  time  domain  difference  (or  recursion)  equation  that  describes  the 
number  sequence  (or  sequence  of  discrete  signal  values)  to  be  modeled  in 
the  z-plane.  In  practice,  the  continuous  functions  in  a  hybrid  system 
are  first  expressed  in  the  s-domain  and  then  transformed  to  the  z-plane 
using  standard  techniques  such  as  partial  fraction  expansion  coupled 
with  table  lookup  or  by  employing  the  inversion  integral  and  contour 


integration.  (The  partial  fraction  expansion  table  lookup  approach  is 
mechanized  in  the  DtSCRET  computer  program.)  On  the  other  hand,  a  dis¬ 
crete  function  (e.g.,  digital  controller)  may  be  first  modeled  with  a 
recursion  equation  and  then  directly  converted  to  the  z-plane  by  substi¬ 
tuting  the  z~n  delay  operator  for  each  discrete  term  in  the  recursion 
equation.  Naturally,  during  the  design  phase,  it  is  the  discrete 
controller  expressed  in  the  z-plane  that  is  first  obtained  and  then  con¬ 
verted  to  a  recursion  equation  using  the  z-n  delay  operator  and  subse¬ 
quently  Implemented  on  a  digital  computer. 

In  analysis  and  design,  no  distinction  is  normally  necessary  between 
the  z-transform  function  derived  from  a  sampled  continuous  element  and 
the  z-transform  function  that  models  a  completely  discrete  element. 
Once  discretized,  all  elements  of  a  hybrid  system  can  be  treated  using 
common  analysis  and  design  techniques  and  tools.  However,  consideration 
must  be  given  to  the  fact  that  discretizing  a  continuous  function  in  the 
z-doraain  only  accounts  for  the  continuous  variables  at  the  sampling 
instance.  In  general,  the  inter-sample  response  Is  not  modeled  with  the 
standard  z-transform.  It  is  necessary  to  investigate  the  inter-sample 
response  using  such  techniques  as  the  continuous  frequency  response  and 
T/N  methods  in  Ref.  1  or  the  advanced  (or  delayed)  z-transform.  Never¬ 
theless,  the  ability  to  model  continuous  and  discrete  elements  in  a 
common  domain  is  one  of  the  most  fundamentally  useful  properties  of  the 
z-transform  (and  tne  w-  or  w'-transf o.*m) . 

B.  FUNDAMENTAL  SAMPLED-DATA  RELATIONSHIPS 

The  fundamental  relationships  for  the  Laplace  transform  of  a  sampled 


signal  are: 

OL 

CT(s)  =  Y.  c(kT)  e_kTs 
k=0 

(1) 
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MUti 


Vtmim'WrrmkL 


CT(s> 


c<p)  t^T^FT  4p 


C'f(s)  -  i  £  c(8-Jf±) 

lf3H»  ■  / 


The  superscript  T  designates  the  time  interval  between  each  sampling 
operation.  These  expressions  are  equivalent  and  each  has  found  varying 
degrees  of  utility  in  sampled-data  or  discrete  system  theory.  All 
assume  that  the  sampling  process  can  be  visualized  as  the  amplitude 
modulation  of  an  impulse  train  6^.(t)  by  the  input  signal  (Fig.  2).  The 
impulse  train  (Eq.  4)  represents  a  series  of  impulses  of  unit  strength 
or  area,  equally  spaced  in  time  and  extending  from  zero  to  plus  Infinity. 


6T(t)  =  6 ( t )  +  5(t  -  T)  +  6(t  -  2T)  + 


-  ]£  fi(t  -  kT) 

k=*0 


The  Laplace  transform  of  $T(t.)  is  given  in  closed  form  as 


sC[5T(t)]  =  1  +  e~sT  +  e“2sT  + 


=  £  e-^T 

k=0 


1  -  e-s^ 


|e“sT|  <  1 


Equation  1  is  the  standard  definition  of  the  z-transform  with  the 
simple  change  of  variable 


z  =  esT 


where  T  is  the  sampling  interval.  Substituting  Eq.  6  into  Eq.  1  con- 
T 

verts  C  (s),  a  nonalgebraic  function  in  s,  to  a  rational  function  in  z. 
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0  T  2T  3T  4T  5T  6T  7T  8T 


figure  2.  Amplitude  Modulation  of  Impulse  Train 
by  Continuous  Signal 
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This  change  of  variable  allows  many  of  the  well-defined  analysis  and 
synthesis  techniques  developed  for  continuous  systems  to  be  applied  more 
readily  to  sampled-data  systems. 

The  relationship  in  Eq.  2  is  obtained  by  exercising  the  method  of 
complex  convolution.  The  dual  relationships  which  are  of  fundamental 
importance  are  stated  below: 

•  The  Laplace  transform  of  the  convolution  of  two 
time  functions  is  equal  to  the  product  of  their 
individual  transforms. 

•  The  Laplace  transform  of  the  multiplication  of 
two  time  functions  is  equal  to  the  convolution  of 
their  transforms  in  the  complex  domain. 

An  analytical  definition  of  the  latter  relationship  is  expressed  as 


,  C  c+j°° 

G(s)  =  2—  I  Gj(p)  G2(s  -  p)  dp 

J  J  c-j=° 


where 


g(t)  =  g^(t)  g2(t) 


oal  <  c  <  Re (s  -  oa2) 


(7) 


(8) 


For  convergence,  the  real  part  of  s  must  be  large  enough  so  that  all  the 
poles  of  -  p)  in  the  p-plane  lie  to  the  right  of  the  poles  of 

Gj(p).  The  abscissa  of  absolute  convergence  of  G|(p)  and  G?(s  -  p)  are, 
respectively,  oaj  and  oa2«  Applying  Eq.  7  to  the  sampling  process 

- ; 1  T *-  1 1  4.  J  —  1  J *-  4  ~  XT  4  1 <■ ~  4  —  - - a-  J - 

lejJlC&euLCU  uy  cue  iauiLipiJ.Gdi.xtui  sjx.  cm  unpuiac  ttdiu  uy  «x  uunLiuuvjua 

time  function  [i„e.,  C  (t)  =  c(t)5^,(t)],  and  recalling  that  the  Laplace 
transform  of  ST(t)  is  given  by  Eq.  5,  results  in  Eq.  2. 


There  are  many  ways  of  deriving  Eq.  3.  Assuming  C(s)  has  at  least 
two  more  poles  than  zeros,  or  the  initial  value  of  c(t)  is  zero  [i.e., 
C(s)  has  a  continuous  impulse  response],  then  the  open  interval  of  inte¬ 
gration  in  Eq.  2  may  be  closed  through  an  infinite  semicircle  in  the 
right  half  plane  as  shown  in  Fig.  3.  The  integral  along  the  infinite 
semicircle  vanishes  as  a  consequence  of  the  assumption  that  the  degree 
of  the  denominator  of  C(s)  is  at  least  two  higher  than  the  numerator 
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Figure  3.  The  Complex  p-Plane 


(Ref.  5).  The  closed  contour  integration  then  reduces  to  the  original 
line  integral  in  F,q.  2.  Cauchy's  integral  formula  then  allows  the 
evaluation  of  Eq.  2  within  the  closed  contour  C2  as  an  infinite  summa¬ 
tion  of  residues  which  include  all  the  poles  of 


1 _ 

1  -  e-T(s-p) 


0 


(9) 


or 

p  =  s  -  ,  k  -  0,  ±1,  ±2,  ...  (10) 


Equation  2  then  reduces  to 


CT(s) 


OO 


-  E 


C(p) 

-  e-T(s-p)] 


(11) 
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P=s-(  j2itk/T) 


The  negative  sign  for  the  summation  is  a  result  of  the  clockwise  contour 
C2*  Evaluating  the  derivative  in  the  denominator  results  in 


4-  fl  -  e~T<>-p)]  =  -Tej2llk  =  -T 

dP  p=s-(j2irk/T) 


(12) 


Equation  11  then  reduces  to 


CT(s ) 


(13) 


If  C(s)  has  a  denominator  one  degree  higher  than  its  numerator  or 
c(0+)  ^  0,  Eq.  13  should  be  modified  to  include  an  additional  initial 
condition  term  (P.ef,  5). 


CT<s)  -  i  ±  c(s  -  Jflt)  +  i  c(0+) 

k=-«  \  / 


(14) 


Restricting  C(s)  to  be  of  order  l/sm,  where  m  >  2  in  Eq.  13  and  m  >  1  in 
Eq.  14,  insures  that  the  infinite  summation  will  be  absolutely  conver¬ 
gent  and  independent  ot  the  order  ot  summation.  However,  the  restric¬ 
tion  on  Eq.  13  can  be  relaxed  to  in  >  1  if  the  sum  is  evaluated  by  taking 
pairs  of  terms  corresponding  to  equal  positive  and  negative  values  of 
the  index  k.  Under  this  condition,  the  sum  in  Eq.  13  will  then  be  abso¬ 
lutely  convergent  (Ref.  11). 

An  alternate  expression  for  Cx(s)  can  be  obtained  from  Eq.  2  by 
closing  the  contour  to  the  left  and  evaluating  the  finite  residues  of 
C(p).  This  contour  avoids  the  problems  of  an  infinite  summation.  For 
this  case,  C(s)  is  required  c ily  to  have  a  denominator  one  degree  higher 
than  its  numerator.  Under  these  conditions,  Eq.  2  reduces  to  the  fol¬ 
lowing  finite  summation  of  residues  corresponding  to  the  poles  of  C(p) 
in  the  p-plane. 
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(15) 


C'r(s)  =  £  residues  ^ 

*■  p=Poles  of  CCp'l 

C.  z-TRANSFORM  AND  THE  INVERSION  INTEGRAL 

The  analytical  derivation  of  a  general  conversion  equation  which 
allows  a  low-rate  discrete  transform  to  be  calculated  from  a  high-rate 
discrete  transform  (TXCONV  computer  program),  relies  on  the  fundamental 
definition  of  the  z-transform  and  application  of  the  discrete  inversion 
integral.  These  two  relationships  are  derived  in  some  detail  in  this 
subsection. 

The  definition  of  the  z-transform  stems  from  the  infinite  summation 

00 

GT(t)  =  E  c(ki')  6(t  -  kT)  ,  k  =  0,  1,  2,  ...  (16) 

k=0 

where  C  (t),  the  sampled  signal,  is  represented  by  the  area  or  strength 
of  impulses  equal  to  the  magnitude  of  c(t)  at  the  sampling  instants 
t  =  kT.  Viewing  the  sampling  process  as  the  amplitude  modulation  of  an 
impulse  train  ST(t)  by  the  signal  c(t)  at  the  sampling  instants  forms 
the  mathematical  basis  for  practical  sampled-data  system  analysis  and 
synthesis.  Such  a  viewpoint  is  justified  if  the  actual  time  during 
which  the  sampler  is  closed  is  short  compared  to  the  time  constants  in 
the  system  under  investigation.  It  is  shown  in  Ref.  5  that  for  a  system 
with  a  single  time  constant  x  =  1/a,  the  error  using  impulse  modulation 
is  less  than  5  percent  for  a  sampler  pulse  width  h  which  is  less  than  or 
equal  to  one-tenth  of  the  time  constant  (i.e.,  h/x  <  1/10).  It  is 
significant  to  note  that  whether  c(t)  is  sampled  physically  or  fic.ti- 
tiously,  or  already  exists  in  pulsed  form,  C  (t)  is  still  representative 
of  an  equivalent  linearized  continuous  signal  c(t)  at  the  sampling 
instants  t  =  kT.  This  point  will  be  elaborated  on  in  the  next  subsec¬ 
tion.  Taking  the  Laplace  transform  of  Eq.  16  produces 


.1 1 


CT(s)  =  c(0)  +  c(T)e_sT  +  c(2T)e~2sT  +  •••  =■  £  c(kT)e-ksT  (17) 

k-0 

In  general,  if  the  Laplace  transform  of  c(t)  is  a  rational  algebraic 
function,  a  closed  form  can  be  found  for  the  infinite  series  represen¬ 
tation  of  C^(s).  The  final  simple  change  in  variable  z  *  es^  results  in 
the  one-sided  z-transform 


CT(z)  *=  £  c(kT)z~k  (18) 

k=0 


For  the  two-sided  z-transform,  the  lower  summation  limit  becomes  minus 
infinity  and  c(t)  is  defined  for  negative  time. 

The  inversion  integral  is  a  closed  form  technique  for  finding  the 
inverse  z-transform  (Eq.  19). 


c(k'r')  -  ~  ^  zk“lcT(z>  dz  (19) 

i  »rt 

Equation  19  is  based  on  the  Laurent  series  expansion  of  F(z)  =  z  Cx(z) 
about  z  =  0.  Expanding  Eq.  18,  the  fundamental  definition  of  the 
z-transform,  produces 

CT(z)  =  c(0)  +  c(T)z~l  +  c(2T)z~2  +  •••  +  c(kT)z“k  +  •••  (20) 

If  we  now  multiply  Eq.  20  by  zk 

F^(z)  =  zk-I.cT(z)  s  c(0)zk-l  +  c(T)zk-2  +  +  c(kT)z~l  +  •••  (21) 

The  desired  output  c(kT)  in  the  Laurent  series  expansion  is  defined  as 
the  residue  of  the  function  FT(z).  This  result  may  be  generalized 
through  application  of  the  Cauchy  theorem  which  states  that  if  the  inte¬ 
gral  FT(z)  is  defined  by 
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(22) 


pT(z)  *  2Sj  fiCi  ^  dz 

and  the  Integral  is  taken  around  a  closed  contour  which  encloses  the 

T 

origin  of  the  z-plane,  then  F  (z)  is  given  by 

Ft(z)  -  0  ,  k  <  -1 

FT(z)  =1  ,  k  **  -1(23) 

FT(z)  -  0  ,  k  >  -1 

where  the  k  =  -1  case  is  recognized  as  the  residue  of  F  (z).  Applying 
this  theorem  term  by  term  to  Eq.  21  results  in  the  discrete  inversion 
integral,  Eq.  ]9.  The  desired  discrete  time  inversion  for  c(kT)  then 
reduces  to  the  practical  evaluation  of  the  residues  of  the  poles  asso- 
dated  with  [z*~  (z))  expressed  in  closed  form  as 


c(kT)  »  residues  of  [z^”J-CT(z)]  at  poles  of  zk~^CT(z)  (24) 
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SECTION  III 


DEVELOPMENT  OF  TRANSFORM  CONVERSION  EXPRESSION 

A.  INTRODUCTION 

The  TXCONV  computer  program  is  mechanized  to  calculate  a  low-rate 
discrete  transform  from  a  given  high-rate  discrete  transform.  This  sec¬ 
tion  analytically  derives  the  transform  conversion  expressions  used  in 
TXCONV.  Detailed  derivations  are  given  in  the  z-,  w-,  and  w'-planes. 
The  mechanization  scheme  in  TXCONV  is  based  on  the  practical  calculation 
of  the  residues  associated  with  a  unique  form  of  the  inversion  integral 
derived  herein.  The  methodology  presented  in  this  section  is  exact  for 
integer  ratios  of  high-to-low  sampling  rate  and  is  based  on  an  equiva¬ 
lent  linearized  continuous  response  for  non-integer  ratios. 

B.  TRANSFORM  CONVERSION  IN  z-DOMAIN 

The  objective  of  the  following  mathematical  development  is  to  derive 
a  closed-form  expression  for  the  low-rate  transform  C^(z)  as  a  function 
of  the  high-rate  transform  CT^(zp)  represented  by 

CT(z)  =  [cT/N(Zp)]T  (25) 

This  transformation  inherently  arises  when  the  output  of  a  system  is 
sampled  at  a  lower  rate  than  the  input  (Fig.  4).  The  superscript  in 
Eq.  25  designates  the  sampling  interval  in  the  z  or  zp  transform.  That 
is. 


CT(z)  - 

—  z  =  esT 

(26) 

CT/N(zp)  - 

-  zp  =  esT/N 

(27) 
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R  y  RT/N 

r  rT/N  rT 

c  y  c  y  <v 
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T/N 

T/N  T 

(phantom) 

Figure  4.  Fast-Input/Slow-Output  Sampling 
with  Phantom  T/N  Output  Sampler 

The  Zp-transform  of  the  sampled  signal  C^N(t)  is  first  calculated  with 
respect  to  a  sampling  interval  of  T/N  producing  the  high-rate  transform 
CT/,N(Zp).  Then,  the  z-transform  of  this  high-rate  transform  is  taken 
with  respect  to  a  T  sampling  interval  producing  the  low-rate  transform 
C  (z).  This  constitutes  the  general  form  of  the  transform  conversion 
addressed  in  this  subsection.  To  simplify  the  notation,  the  z  and  zp 
designation  will  he  suppressed  and  C~  and  will  be  used  to  designate 

the  z  and  zp  transforms. 

We  will  assume  that  the  sampling  ratio  N  in  Eq.  25  can  be  any  inte¬ 
ger  or  non-integer  rational  value.  However,  only  integer  values  of  N 
are  allowed  in  most  practical  sampled-data  systems.  More  will  be  said 
about  this  in  the  next  subsections.  For  the  present,  we  proceed  with 
the  derivation  for  rational  values  of  the  sampling  ratio  and  subse¬ 
quently  treat  integer  values  as  a  special  case  of  the  more  general  non¬ 
integer  case. 

The  respective  sampled  signals  in  Eq.  25  are  defined  in  the  z-domain 
using  Eq.  18. 

CT/N  =  £  c(kT/N)z"k  ,  zp  =  esT/N  (28) 

k=0 

oo 

CT  =  £  c(kT)z“k  ,  z  =  esT  (29) 

k=0 
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Transforming  Eq.  28  back  Into  the  tLme  domain  using  the  inversion  inte¬ 
gral  (Eq.  19)  results  in 


c(kT/N) 


(30) 


It  is  important  to  recognize  that  although  Eq.  30  provides  information 
only  at  discrete  instances  of  time  separated  by  T/N  seconds,  a  linear¬ 
ized  continuous  time  function  c(t)  can  be  obtained  from  the  solution  of 
Eq.  30  by  replacing  kT/N  with  t.  This  linearized  system  response  agrees 
with  the  sampled  response  at  the  sampling  instants  t  =  kT/N.  Moreover, 
this  linearized  response  also  exactly  characterizes  the  low-rate  sampled 
response  c(kT)  with  t  replaced  by  kT  for  integer  values  of  N.  It 
approximates  the  low-rate  sampled  response  c(kT)  for  non-integer  values 
of  N  by  assuming  that  c(kT/N)  is  the  high-rate  sampled  response  of  a 
continuous  system  c(t).  For  example,  if  C(s)  contains  only  simple  poles 
at  (aj,  a2»  a-j,  ..,),  the  general  closed-form  discrete  time  solution 
from  Eq.  30  is  represented  by 

c(kT/N)  =  Aie"aikT/N  +  A2e“a2kT/N  +  A3e~a3kT/N  +  •••  (31) 

where  (A,  A£,  A3,  . ..)  represent  the  residues  of  the  simple  poles.  If 
kT/N  is  replaced  by  t  in  Eq.  31,  the  linearized  continuous  response  is 
obtained  (Eq,  32)  which  exactly  matches  the  sampled  response  c(kT/N)  at 
t  *  kT/N. 


c(t)  “  A3e  +  A2e  a2t  +  A3e  a3t  +  (32) 

It  is  readily  apparent  that  sampled  responses  at  other  than  the  T/N 
interval  can  be  obtained  from  Eq.  30  via  a  change  in  the  ordering 
index  k.  For  the  T  sampling  interval,  substituting  k  =  kN  in  Eq.  30 
produces 
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c(kT) 


pT/N  kN-l  , 

C  zp  dzp 


Substituting  Eq.  33  into  Eq.  29  results  in  Eq.  34. 


CT  = 


y  _L  CT/N  kN-l  -k 

£o  2"J  p  PJ 


Interchanging  the  summation  and  integration  in  Eq.  34  produces  Eq.  35. 


C?  - 


^  d>  cp«  e  (A-!)*  l 

J  -m:,  k=n  p  2p 


The  infinite  summation  in  Eq.  35  is  recognized  as  a  geometric  progres¬ 
sion  in  (zNz~ll  which  can  be  placed  in  closed  form  as  indicated  in 
P 

Eq.  36. 


cT  • 


1  —  z^z  ^ 


Equation  36  is  the  final  result  which  can  be  used  to  calculate  any 
general  low-rate  z-transf orm  from  a  given  high-rate  Zp-transf orm.  To 
evaluate  Eq.  36,  the  integration  contour  in  the  Zp-plane  can  be 
selected  to  include  all  the  poles  of  c^^/z  and  exclude  the  N  poles  of 
(zNz'l).  Alternatively,  the  contour  can  include  only  the  N  poles  of 
(z^z^l)  and  exclude  the  poles  associated  with  C^^/’.p.  For  either  ap¬ 
proach,  the  problem  reduces  to  calculating  the.  residues  of  the  enclosed 
poles.  The  computationally  more  convenient  method  is  to  evaluate  the 
residues  of  CT/N/zp.  Equation  36  then  reduces  to  the  finite  summation 
given  in  Eq.  37. 
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residues 


Zp  2  “  2^ 

P  Zp*Poles  of  CT/N/zp 


where 


CT  -  CT(z) 


Low-rate  transform,  z  =  esT 


CT/N  .  CT/N(zp)  High  -rate  transform,  Zp  =  esT/N 

C.  TRANSFORM  CONVERSION  IN  w-  AND  w '-PLANE 

The  bilinear  transformation  from  the  z-plane  to  the  w-  and  w'-planes 
are  defined  by; 


1  +  w 
I  -  w 


z  -  I 
z  +  1 


2/T  +  w' 
2/T  -  w' 


2  z  -  1 
T  z  +  I 


Since  w'  is  related  to  w  by  a  scale  factor  2/T,  the  bilinear  transforma¬ 
tion  can  be  expressed  as 


A  +  w 
A  -  w 


where  w  represents  either  w  or  w'  and  A  =  1  for  w  and  A  =  2/T  for  w'. 
Substituting  Eq.  40  into  Eqs.  28  and  29  produces 


w  4  4-  w 

E  ««/»)  £_13l 

k=0  P  PJ 


.  zp  “  1 
WP  AP  zp  4-  1 


CT  =  f;  c(kT)  ' 

k=0  A  w. 


z  -  1 
z  +  1 
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where  A  is  associated  with  the  low-rate  transform  and  Ap  with  the  high- 
rate  transform.  Transforming  Eq.  30  into  the  w-  or  w'-plane  involves  a 
change  in  the  integrtion  variable  given  by 


dZp  =  (Ap  -  wp)2  dWp 


Substituting  Eqs.  40  and  43  into  Eq.  30  and  simplifying  produces 


c(kT/N)  = 


An  +  Wr 


CT/N  "P  _  P 


Ap  -  Wp  Ap  -  Wp  Ap  +  Wp 


Changing  the  sampling  index  to  k  =  kM  in  Eq.  44,  substituting  Eq.  44 
into  Eq.  42,  and  interchanging  the  summation  and  integration  results  in 


y  /ap  +  wp\N  M  +  w\~1  k  \  2Ap  dwp 

k~0  Up  -  wp /  U  -  w/  |  Ap  -  Wp  Ap  +  Wp 


Placing  the  infinite  summation  in  closed  form  reduces  Eq.  45  to 


ct  =  .J__  (f>  ct/n _ i _  . .  2a£ _  _  d..ri, 

2nj  )rcl  1  -  X  Ap  -  wp  Ap  +  wp 


where 


/Ap  +  wp\N  /A  +  w\- 

Up  "  wp/  \  A  -  w  / 


I  X|  <  1 
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As  was  done  in  Eq.  36,  a  closed  contour  Cj  is  taken  which  includes  all 
the  poles  of  the  integrand  in  Eq.  46  except  the  poles  associated  with 
the  infinite  summation  term  1/1  -  X.  The  free  1/Ap  -  wp  terra  must  also 
be  excluded  since  it  is  in  the  same  contour  region  as  the  1/1  -  X  term. 
It  is  required  that  the  infinite  sum  in  Eq.  45  be  absolutely  convergent 
(i.e.,  satisfy  Eq.  48)  over  the  region  of  integration.  For  then,  the 
integral  remains  finite  and  the  summation  of  integrals  in  Eq.  45  equals 
the  integral  of  the  summation.  This  condition  assures  a  closed  form 
solution  for  Eq.  46.  The  integral  then  reduces  to  a  summation  of  resi¬ 
dues  given  by: 


CT  = 


where 


Wp=Po]es  of  CT'/N/(Ap+Wp ) 

(49) 

Wp  =  High-rate  transform  variable 
w  *  Low-rate  transform  variable 


y,  residues 


Ap  +  Wp  1  X  Ap  Wp 


w  *  w',  Wp  =  Wp  for  Ap  =  2N/T  and  A  =  2/T 

w  *  w,  Wp  =  Wp  for  Ap  =  1  and  A  =  1 


D.  GENERAL  MULTI-RATE  CONFIGURATIONS 


In  the  preceding  derivation  leading  to  Eq.  36  (and  Eq.  49)  it  was 
shown  Lli-.L  the  high  to  low-rate  sampling  ratio  N  can  be  any  rational 
value.  There  are  three  general  multi-rate  cases  of  particular  interest 
in  sarapled-data  systems: 


®  If  N  is  an  integer,  the  result  obtained  from 
F.qs.  36  or  49  are  exact  for  all  multi-rate  con¬ 
figurations,  This  is  easily  seen  since  the.  exact 
discrete  low-rate  information  c.(kT)  can  be  selec¬ 
tively  extracted  from  the  original  high-rate  dis¬ 
crete  signal  c(kT/N).  Therefore,  no  use  is  made 
of  the  linearized  continuous  response  function. 
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•  For  non-integer  values  of  N,  the  results  from 
Eqs.  36  or  49  are  exact  if  the  original  high-rate 
discrete  response  c(kT/N)  is  the  sampled  response 
from  a  completely  continuous  linear  function  c(t) 
or  C(s). 

•  If  N  is  a  non-integer  and  c(kT/N)  is  the  response 
from  a  high-rate  discrete  function  (e.g.,  digital 
computer  algorithm),  the  results  from  Eqs.  36  or 
49  are  at  best  approximate  since  the  low-rate 
response  c(kT)  is  based  on  a  linearized  continu¬ 
ous  model  of  the  original  discrete  function. 

The  most  practical  multi-rate  configurations  are  those  which  can  be 
analyzed  with  integer  values  of  N.  The  high-rate  transform  can 

originate  from  a  completely  continuous  function,  a  completely  discrete 
function,  or  any  combination  of  continuous  and  discrete  functions,  and 
Eqs.  36  or  49  provide  the  exact  low-rate  discrete  transform  for  integer 
ratios  of  sampling  rate. 

Consider  the  fundamental  fast-input/slow-output  multi-rate  sampling 
configuration  in  Fig.  5  where  M  represents  the  transfer  function  of  a 
data  hold  and  G  a  continuous  system  in  the  s~plane„  In  general,  the 
output  sampling  interval  T2  is  greater  than  the  input  sampling  interval 
Tjj  however,  a  direct  integer  ratio  between  output-to-input  sampling  is 
not  necessarily  implied.  This  simple  multi-rate  configuration  repre¬ 
sents  a  complex  (but  practical)  situation  which  can  be  analyzed  using  a 
common  sampling  period  T  such  that 


T 

M  "  T1 

T 

N 

=  *2 

M,N  =  integer 

(50) 

For  example,  if 

the 

input 

is  sampled  at  20  cps 

and  the  output  at 

13.33  cps;  Tj  =  0. 

050, 

t2  =  0 

.075,  and  a  choice  of  T 

*  0.150  produces 

T 

3  ~  Ti 

T 

2 

-  T2 

T  =  0.150 

(51) 

The  multi-rate  configuration  in  Fig.  5  then  reduces  to  the  general  form 
in  Fig.  6.  The  output  equation  foi  lig.  6  is  given  by 
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Figure  5.  Fundamental  Fast-Input/Slow-Output 
Multi-Rate  Sampling,  T2  >  T^ 


Figure  6.  Fast-Input/Slow-Output  Sampling 
with  Common  Sampling  Period  T 


CT/N  „  [gmrT/M]T/N 


(52) 


Computationally,  Eq.  52  is  more  involved  than  a  general  slow-input/fast- 
output  system,  since  the  T/N  operator  does  not  "operate  through"  the 
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interval  is  greater  than  the  inner  sampling  interval  T/M.  Moreover,  T/N 
is  not  necessarily  an  integer  multiple  of  T/M,  which  further  complicates 
the  analysis.  Fortunately,  we  are  free  to  add  a  mathematical  or  phantom 
sampler  to  the  output  (Fig.  7)  which  operates  at  an  integer  multiple  of 
the  output  sampling  rate  or  at  a  submultiple  of  the  output  sampling 
interval  T/N.  This  mathematical  convenience  is  valid  since  the  actual 
output  sampler  T/N  simply  rejects  all  the  unwanted  samples  from  the 

phantom  sampler  T  .  This  cimplif ication  overcomes  the  above  complica- 

* 

tions  and  facilitates  a  solution  to  Eq.  52  using  the  high-rate  to  low- 
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Figure  7.  Fast-Input/Slow-Output  Sampling  with 
Phantom  Sampler 


rate  transform  expression  in  Eqs.  36  or  <^9.  With  the  additional  sampler 

T  ,  the  output  equation  from  Fig.  7  becomes 
* 


CT/N  =  [(GM)T*  ^T/MjT/N  (53) 

The  problem  is  now  reduced  to  finding  the  individual  transforms  (GM)T* 

and  Rt/m  (and  their  product)  using  a  common  definition  for  the  transform 

variable  z,  w,  or  w'.  For  the  more  general  case  of  where  M/N  is  not  an 

integer,  a  value  for  T  must  be  selected  such  that  T  is  both  smaller 

*  * 

than  and  an  integer  submultiple  of  T/M  and  T/N.  An  obvious  choice  is 

T  =  T/MN;  however,  the  largest  compatible  value  for  T  is  composed  of 
*  * 

the  prime  factors  of  the  integer  product  MN.  The  smaller  T  is,  the 

* 

higher  the  order  of  the  numerator  and  denominator  polynomials  in  the 
rT/m  term.  This  increase  in  order  is  the  result  of  substituting  the 
higher  rate  transform  variable  (e.g.,  z  •-  esT*)  associated  with  the 

j|f 

(GM)r*  term  into  the  lower  rate  R^M  term  (z0  =  esT^)  to  form  an  over¬ 
all  high-rate  discrete  transform  product.  For  example,  if  z  = 

sT  /  3  * 

and  zm  =  e  '  ,  then  the  zm  transform  variable  can  be  defined  by 

zm  =  z2  and  the.  (GM)^^R^^  product  can  be  formed  using  the  common 

transform  variable  z  .  Therefore,  a  T  composed  of  the  prime  factors  of 

*  * 

MN  reduces  to  a  minimum  the  resulting  order  of  the  term.  This  in 

turn  may  reduce  the  computations  required  to  calculate  the  low-rate  T/N 
transform  in  Eq.  53  (using  Eqs.  36  or  49)  since  the  number  of  residues 
in  the  (GM)l*RT/M  product  will  be  at  a  minimum.  Using  the  prime  nota¬ 
tion,  the  general  sampling  structure  in  Eq.  53  becomes 
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CT/N  =  [(GM)1/^*  rT/M]T/N 


(54) 


Inserting  the  previous  numerical  sampling  values  into  Eq.  54  will 
help  summarize  and  clarify  the  general  results. 


cT/2  =  [(Gm)T/6  rT/31t/2  (55) 

The  selection  of  =  T/6  allows  the  formation  of  the  inner  transform 
product  (GM)Ty/kRT/3  an(j  at  the  same  time  provides  an  integer  ratio  of 
outer-to-inner  sampling  periods.  Equation  55  is  now  in  a  form  that  can 
be  easily  solved  by  Eqs.  36  or  49. 

A  second  fundamentally  important  system  configuration  is  shown  in 
Fig.  8.  Here,  G^'/N  represents  a  discrete  model  of  the  digital  computa¬ 
tions  in  a  computer  (e.g.,  digital  control  laws).  is  a  data  hold 

device  that  models  the  holding  of  information  in  a  storage  register 
between  sampling  intervals.  The  actual  computational  time  in  the  com¬ 
puter  is  assumed  negligible  in  this  case  and  is  not  considered.  How¬ 
ever,  computational  delays  can  be  easily  handled  with  appropriate  delay 
factors  and  the  advanced  z~,  w-,  or  w '-transforms.  Simple  algebraic 
signal  flow  tracing  produces  Eq.  56,  the  output  equation  for  Fig.  8. 

CT  =  [MnGT'/NRT/'NjT  (56) 


Figure  8.  Fast-Input/Slow-Output  Sampling  with 
Discrete  System  Component 
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The  same  situation  persists  here  as  in  Eq.  52  and  we  are  unable  to 
operate  through  with  T.  Adding  a  phantom  sampler  to  the  output  allows 
us  to  write 


CT  =  [MT/NGT/NRT/N]T  <57) 

The  CT  transform  in  the  2-plane  can  now  be  easily  obtained  from  Eq.  36 
by  evaluating  the  residues  of  the  inner  transform  product  with  the 
transform  variable  defined  as  z  =  esT//^. 
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SECTION  IV 


DESCRIPTION  OF  DISCRET  COMPUTER  PROGRAM 


A.  INTRODUCTION 

DISCRET  is  a  versatile  and  relatively  accurate  digital  computer  pro¬ 
gram  that  transforms  a  continuous  s-plane  transfer  function  into  a  valid 
discrete  transfer  function  in  the  s-,  w-,  or  w'-plane.  The  program  cal¬ 
culates  the  discrete  transformation  for  a  system  of  the  general  form 
depicted  in  Fig.  9. 


Figure  9.  General  Form  of  Sampled  System 


M(s)  in  Fig.  9  represents  the  s-plane  transfer  function  model  of  a 
data  hold,  G(s)  the  s-plane  representation  of  a  continuous  system,  and 
eaiS  the  time  factor  in  the  advanced  or  delayed  discrete  transform.  For 
the  standard  z-,  w-,  or  w'-transf orm,  AT  ■  0,  The  output  equation  for 
Fig.  9  is 

CT  =  [eATs  G(s)M(s)]^  rT  (58) 

The  superscript  T  denotes  the  sampling  period  in  the  z~,  w-,  or  w'- 
transform  (i.e.,  z  =  es^).  The  DISCRET  program  calculates  the  general 
transform  given  by 


26 


(59) 


[ e^Ts  G(s)M(s)]a 


This  general  transform  is  present  in  most  open-loop  and  closed-loop 
sampled-data  or  discrete  systems. 

DISCRET  is  written  in  FORTRAN  for  the  Control  Data  Corporation  (CDC) 
CYBER  175  series  computer.  It  can  handle  pole  multiplicity  up  to  three 
and  system  order  up  to  50th.  Double  precision  arithmetic  is  used 
throughout  the  program. 

B.  PROGRAM  OPTIONS 

DISCRET  calculates  the  practical  discrete  transformations  required 
to  analyze  and  design  realistic  sampled-data  systems.  The  program  can 
execute  any  combination  of  the  following  options: 

•  Transform  Options 

-  z-transform 

-  w-transform 
w' transform 

•  Data  Hold  Options 

-  None 

-  Zero  order  hold 

*11,..*.  _-.J_  —  L  _  1  J 

“  rubl  uiuct  uuxu 

-  Second  order  hold 
Slewer 

•  Time  Increment  Option 

-  Standard  transform 
Delayed  transform 
Advanced  transform 
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C.  MATHEMATICAL  DESCRIPTION  OF  PROGRAM 


The  computer  code  in  DISCRET  is  based  on  the  computer  program  in 
Ref.  12.  Extensions  and  modifications  have  been  made  to  accommodate 
additional  options.  The  analytical  basis  for  the  computer  algorithms  is 
discussed  in  Ref.  12.  However,  to  provide  a  complete  description  of  the 
program,  some  of  the  details  are  repeated. 

Consider  the  partial  fraction  expansion  of  the  general  expression  in 
Eq.  59. 


[eATsG(s)M(s) ]T 


Ft(z) 


_ *1  +  *21  +  *22 

(s  +  ai)  (s  +  a2)  (s  +  a2)2 


+ 


*31 


*32 


*33 


(s  +  a3)  (s  +  a3)z  (s  +  a3)- 


(60) 

Ft(z)  in  Eq.  60  is  a  z-plane  term  determined  by  the  data  hold  se¬ 
lected.  For  example,  for  a  zero  order  hold  (ZOH),  Eq.  59  becomes 


f  eATsG(s)M(s) 1T 

L  *  '  '  J 


LzJ—  *atsg(s-) 


Z  -  1 


,ATs 


G(s)  T 

s  J 


(f>  i  i 


where  the  ZOH  introduces  a  pole  at  8  =  0  and  a  z-plane  factor  z  -  1/z. 
For  other  data  holds,  a  similar  situation  exists.  This  can  be  readily 
verified  by  applying  the  data  hold  transfer  functions  in  Table  1  to 
Eq.  59. 
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28 


nit  rill  tute. 


TABLE  1.  DATA  HOLDS  IMPLEMENTED  IN  DISCRET 


DATA  HOLD 

TRANSFER  FUNCTION 

Zero-Order 

Hold 

1  _  e-sT 

M0  -  :  • 

First-Order 

Hold 

"1  ■  Mo(s  +i) 

Second-Order 

Hold 

h2  -  »s(=2  +Ir  s  +  72) 

Slewer 

Hold 

«sU„  -  ^ 

In  DISCRET,  the  individual  partial  fraction  expansion  terms  in 
Eq.  60  are  first  transformed  into  the  w-plane  using  the  following 
advanced  w-transforras : 


eATs 
s  +  a 


T 


e~aAT  |  _ I  +  w _ ) 

1  +  e-a^  j  w  +  (l  -  e-a^/l  +  e“aT)  j 


(62) 


,ATs 


(s  +  a)2 


ATe~aAT  l _ 1  +  w _ | 

1  +  e_a^  w  +  (l  -  e~a^/l  +  e~aTj  j 


Te~aTe~aAT  ) _ (1  +  w)(l  -  w) _ } 

( 1  +  e“aTJ^  (  [w  +  (l  -  e“a^/l  +  e-a^} ] ^  ! 
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r  eATs  ]T 

(AT)2e_aAT 

1  1  +  w  ) 

(6  +  a)^J 

2(l  +  e-aT) 

[  w  +  (l  -  e~a^/l  +  e-aT)  ) 

(1+  w)( 1  -  w) 


+  T2(l  +  2A)e-aTe~-aAT  ( _ 

2(l  +  e~a^J^  (  [w  +  (l  -  e-a^/l  +  e~aT)]^ 


T2e-2aTe-aAT  j _ (1  ■+  w)(l  -  w)2 _ | 

(it-  e~aT}^  (  [w  +  ( 1  -  e~aT/l  +  e'aT)  ]  ^  j 


(64) 


In  Eqs.  62-64, 

T  =  Sampling  interval  (sec) 

AT  =»  Time  advance  (sec),  0<AT<T,  0<A<1 

The  w-plane  is  related  to  the  z-plane  by  the  bilinear  transformation 
in  Eq.  65. 


W  ~ 


Z  ~  1 
Z  +  1 


1  +  w 
1  -  w 


2  »  ^sT 


(65) 


The  corresponding  z-plane  transforms  for  the  partial  fraction  expansion 
terms  in  Eq.  60  are  shown  in  Eqs.  66-68  (Ref.  2). 


r  *  m _  IT 
eui& 

s  +  a 


-nAT 

e 

z  -  e~a"^ 


»ATs 


(s  +  a)2 


ATe-a^^z  Te~aTe~a^Tz 


z  —  f 

*  e  [z  -  e 


-  P-aT)2 


(67) 


,ATs 


(s  +  ap 


(AT)2e-aATz  T2(l  +  2A)e-aTe“afiTz  T2e~2aTe“aATz 

*+* - -  11  — —  ■  -f- 


2(z  -  e~aT)  2[z  -  e-aT)^  (z  -  e“aT)' 


(68) 
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The  s-plane  expressions  in  Eqs.  62-64  and  66-68  represent  Laplace 
transforms  of  functions  advanced  in  time.  For  example,  e^s/(s  +  a)  is 
the  Laplace  transform  of  the  continuous  time  function  e-at  advanced  by 
AT  seconds.  That  is, 

e-a(t+AT)  (  0  <  (t+AT)  <  (t+T)  (69) 


Sampling  the  advanced  time  function  in  Eq.  69  with  period  T  and  taking 
the  2-transform  results  in  the  advanced  z-transform 


(70) 


and  the  advanced  w-transform 

e-aAT  (  _ w  -t  1 _ | 

1  +  e-aT  |  w  +  (l  -  e_aT/l  +  e”aT)  j 


(71) 


Numerical  calculations  in  DISCRET  are  carried  out  in  the  w-plane  to 
improve  the  accuracy  of  the  cross-multiplications  necessary  to  form  the 
numerator  of  the  discrete  transfer  function.  The  discrete  numerator  is 
formed  by  multiplying  each  partial  fraction  expansion  numerator  by  all 
denominator  terms  except  its  own  and  then  summing  the  resultant  prod¬ 
ucts.  In  general,  the  poles  of  the  z-transfer  function  tend  to  migrate 
towards  the  unit  circle  in  the  z-plane  (i.e.,  z  ■*>  1).  Computationally, 
severe  loss  of  accuracy  can  result  from  the  summation  of  individual  par¬ 
tial  fraction  expansion  terms  in  the  z-plane.  This  inherent  inaccuracy 
can  be  minimized  by  performing  all  possible  calculations  in  the  w~  or 
w'-plane  where  the  poles  are  more  reasonably  separated.  Therefore,  the 
w-plane  Eqs.  62-64  are  implemented  in  the  computer  code.  It  then 
becomes  a  simple  task  to  calculate  the  corresponding  z-  and  w'-trans- 
forras  using  the  bilinear  transformations  in  Eq.  72. 


The  s-plane  partial  traction  expansion  terms  in  Eq.  60  are  obtained 
in  the  subroutine  PARTFR.  These  are  passed  to  the  subroutine  WPLN  which 
calculates  the  w-piane  transforms  via  Eqs.  62-64.  The  time  factor  term 
e  represents  a  time  advance  of  less  than  one  sampling  interval  T 
(i.e.,  0  <  AT  <  T).  For  time  delays,  an  additional  delay  factor 
z’1  =  1  -  w/l  +  w  is  added  to  Eqs.  62-64  and  AT  is  defined  by 

AT  =  1  -  (D/T) 

T  =  Sampling  interval  (sec) 

D  *  Delay  (sec 

In  this  manner,  the  program  can  calculate  either  the  advanced  or  delayed 
discrete  transform.  The  user  inputs  a  positive  time  advance  (AT),  a 
negative  time  delay  (D) ,  or  zero  for  the  time  increment.  From  this 
information  alone,  the  program  calculates  the  advanced,  delayed,  cr 
standard  discrete  transform. 

The  computer  code  automatically  inserts  a  user-selected  data  hold 
(Table  1)  Into  the  implementation  scheme.  The  appropriate  s-plane  zeros 
and  poles  associated  with  each  data  hold  is  inserted  into  the  s-plane 
continuous  system  G(s)  (Eq.  60)  by  the  main  program  module  ADVANZ.  The 
z-plane  term  F  (z)  in  Eq.  60  is  added  to  the  transformation  by  the  sub¬ 
routine  WPLN. 

D.  PROGRAM  STRUCTURE 

The  basic  structure  of  DISCRET  consists  of  the  ADVANZ  main  program 
and  two  primary  subroutines  PARTFR  and  WPLN.  These  three  program 
modules  along  with  their  supporting  subroutines  are  shown  in  Fig.  10. 
The  main  program  ADVANZ  first  calls  PARTFR  to  obtain  the  partial  frac¬ 
tion  expansion  terms  in  the  s-plane  and  then  calls  WPLN  to  execute  the 
conversion  to  the  z-,  w-,  or  w'-plane.  No  external  libraries  are  used 
with  the  exception  of  the  normal  system  routines  that  support  FORTRAN. 
Parameters  are  passed  between  the  main  program  and  the  two  primary 


3  2 


subroutines  entirely  in  a  COMMON  data  structure.  This  lack  of  formal 
parameter  passing  via  subroutine  arguments  was  initiated  to  allow  the 
program  to  be  easily  converted  to  an  overlay  structure  in  the  TOTAL 
(Ref.  13)  computer  program  at  AFWAL/FIGC.  This  permits  interactive 
operation  of  DISCRET  as  an  option  in  TOTAL.  The  overlay  version  at 
AFFDL  transfers  the  main  program  functions  to  the  TOTAL  main  overlay. 
The  main  program  ADVANZ  is  then  treated  as  a  primary  overlay  with  sub¬ 
routines  PARTFR  and  WPLN  converted  to  secondary  overlays. 

The  source  listing  for  DISCRET  in  Volume  III  is  set  up  in  a  standard 
program-subroutine  structure  (i.e.,  a  main  program  followed  by  its  sub¬ 
routines).  This  listing  also  includes  (in  the  comment  code)  the  re¬ 
quired  changes  to  run  the  program  in  an  overlay  structure.  Dividing  the 
program  into  overlays  reduces  the  amount  of  computer  memory  required  to 
execute  the  program.  The  overlay  code  Is  highlighted  with  a  star  (*) 
character  is  column  one.  Removing  this  code  from  comment  will  allow 
overlay  operation.  To  complete  this  turnover,  the  main  program  card  for 
ADVANZ  and  the  subroutine  cards  for  PARTFR  and  WPLN  must  be  deleted.  In 
addition,  the  two  call  subroutine  statements  located  in  ADVANZ  must  also 
be  deleted. 

DISCRET  can  be  run  in  either  a  batch  or  interactive  mode.  The 
source  code  in  Volume  HI  is  set  up  for  batch  operation.  The  prompting 
code  to  run  the  program  in  an  interactive  mode.  Is  also  Included  in  the 
comment  section  of  the  main  program  ADVANZ.  This  code  can  be  identified 
by  the  charac  ters  "CINT"  in  the  first  four  columns. 

E.  DESCRIPTION  OF  SUBROUTINES 

This  subsection  presents  a  brief  description  of  the  routines  used  in 
DISCRET.  The  general  program  structure  contains  a  main  program  ADVANZ, 
two  primary  subroutines  PARTFR  and  WPLN,  and  16  supporting  subrou¬ 
tines,  These  routines  are  all  coded  in  FORTRAN. 


34 


1 •  Program  ADVAN2 


a.  Purpose 

This  is  the  main  program  for  DISCRET.  It  reads  the  input  from  data 
cards  and  transfers  it  to  internal  variables  and  arrays  that  are  located 
in  a  common  data  structure.  The  program  adds  appropriate  poles  and 
zeros  to  the  input  s-plane  transfer  function  according  to  the  data  hold 
selected.  It  calls  PARTFR  to  calculate  the  partial  fraction  expansion 
terms  and  then  calls  WPLN  to  execute  the  discrete  transformations  to  the 
z-,  w-,  or  w' -plane. 

b.  Input/Output 

All  data  are  read  from  data  cards  and  transferred  to  PARTFR  and  WPLN 
via  labeled  common. 

2.  Subroutine  PARTFR 


a.  Purpose 

This  routine  takes  the  plant  description  in  terras  of  poles  and  zeros 
and  outputs  the  partial  fraction  expansion  coefficients  corresponding  to 
each  pole.  The  program  drops  identical  poles  and  zeros  and  then  calcu¬ 
lates  the  polynomials  needed  for  evaluating  the  partial  fraction  expan¬ 
sion  terms. 

b.  Input/Output 

All  inputs  and  outputs  for  this  subroutine  are  handled  entirely  by 
common  statements. 

3 .  Subroutine  WPLN 

a.  Purpose 

This  subroutine  takes  the  partial  fraction  expansion  coefficients 
and  the  corresponding  poles  and  calculates  the  w-plane  transformation. 
After  determining  the  w-plane  numerator  and  denominator  for  each  partial 
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fraction  expansion  term,  the  cross  product  is  formed.  Each  numerator  is 
multiplied  by  all  denominators  except  its  own.  The  resultant  polynomi¬ 
als  in  w  are  then  summed  by  adding  coefficients  for  each  power  of  w. 
The  subroutine  then  calls  a  root  solver  to  find  the  zeros  of  the  numera¬ 
tor  of  the  overall  transfer  function.  The  roots,  poles  and  w-plane 
polynomials  may  then  be  transformed  to  the  z-  or  w'-plane, 

b.  Input/Output 

All  inputs  and  outputs  for  this  subroutine  are  handled  entirely  by 
common  statements. 

4.  Subroutine  CONVRT  (A,  NA,  PN,  NPN) 


a.  Purpose 


The  purpose  of  this  routine  is  to  change  the  format  of  an  array  with 
real  polynomial  coefficients  and  corresponding  powers  of  s  to  allow  a 
place  for  the  imaginary  part  of  the  coefficient. 


b.  Input 

1)  A:  A  double  precision  arrav  with  a  two  place 

format  such  that  each  real  coefficient  of  a  poly¬ 
nomial  is  immediately  followed  with  its  corre¬ 
sponding  power  of  s. 

2)  NA:  Number  of  occupied  locations  in  array  A. 


c.  Output 

1)  PN:  A  double  precision  array  with  a  three  place 
format  such  that  each  real  coefficient  is  fol¬ 
lowed  by  a  zero  in  the  next  location  and  the 
corresponding  power  of  s  in  the  third  location. 

2)  NPN:  The  number  of  occupied  locations  that  are 

used  in  PN(3*NA/2). 
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5.  Subroutine  DIRIV  (PN,  NPN,  PD,  NPD) 


a.  Purpose 

This  routine  takes  the  derivative  of  a  transfer  function  with  the 
numerator  polynomial  located  in  PN  and  the  denominator  polynomial 
located  in  PD.  It  then  stores  the  numerator  of  the  derivative  in  PN  and 
the  denominator  in  PD.  If 


then 


and 


P  =  V]  Aisn-i 
i=0 


j-  (P)  =  £  <n  -  i)AiSn-i-l 

dS  1=0 


PN  “  ( PD)-Fbj 

PD  =  PD*PD 


b.  Input 

1)  PN:  A  double  precision  array  containing  the 

numerator  polynomial  in  the  format:  real  part, 

imaginary  part,  and  the  order  of  s  stored  in 
back-to-back  locations. 

2)  NPN:  Number  of  occupied  elements  in  array  PN. 

3)  PD:  A  double  precision  array  containing  the. 

denominator  polynomial  in  the  same  format  as  the 
numerator  polynomial. 

4)  NPD:  Number  of  occupied  elements  in  array  PD. 
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c.  Output 


1)  PN:  Numerator  polynomial  (in  the  same  format  as 
the  input)  of  the  derivative  function. 

2)  NPN:  Number  of  occupied  elements  in  PN. 

3)  PD:  Denominator  polynomial  of  the  derivative 

function. 

4)  NPD;  Number  of  occupied  elements  in  PD. 

6.  Subroutine  MLTPL  (C,  N,  D,  E,  M) 

a.  Purpose 

This  routine  multiplies  the  real  polynomial  coefficients  by  a  scale 

factor  and  stores  the  resultant  polynomial  in  a  new  array. 

b.  Input 

1)  C:  A  Double  precision  array  containing  poly¬ 

nomial  coefficients  and  corresponding  powers  of 


2)  N:  Number  of  occupied  elements  in  array  C. 

3)  D:  Scale  factor  which  multiplies  all  odd  loca¬ 

tions  of  array  C. 


c.  Output 


1  \  T7.  A  J 1.  T  ^  J  r.  V ~ 1  ~  J 

i  /  iu  •  n  uuuo-tc  ptct-iaiuu  at  tay  ocaicu 

polynomial  coefficients. 

2)  M1;  Number  of  occupied  elements  in  array  E. 


7.  Subroutine  FORK  (RD,  NRD.  MULT,  NEGLCT,  P2,  NP2 ) 


a.  Purpose 


The.  purpose  of  this  routine  is  to  form  the  denominator  polynomial 
that  will  be  used  to  evaluate  the  partial  fraction  expansion  coefficient 
corresponding  to  a  pole  of  the  plcnl.  The  routine  multiplies  all  of  the 


poles  together,  excluding  the  pole  (and  its  conjugate  if  the  pole  has  an 
imaginary  part)  for  which  the  partial  fraction  expansion  coefficient  is 
being  sought. 

b.  Input 

1)  RD:  A  double  precision  array  containing  the 

poles  of  the  plant,  real  part  then  imaginary 
part. 

2)  NRD:M)S:  Two  times  the  number  of  distinct  poles 
of  the  plant  (number  of  locations  used  in  array 
RD). 

3)  MULT:  Array  contains  multiplicity  corresponding 
to  each  pole  contained  in  array  RD  (NRD/2  loca¬ 
tions  ) . 

4)  NEGLCT:  An  integer  array  containing  all  zeros 

except  in  the  location  corresponding  to  the  pole 
for  which  the  partial  fraction  expansion  coeffi¬ 
cient  is  being  determined,  where  a  one  appears. 

If  the  pole  is  complex,  a  one  also  appears  in 
location  corresponding  to  the  conjugate  of  the. 
pole  (NRD/2  locations). 

c.  Output 

1)  P2;Utfb:  A  double  precision  array  containing  a 

polynomial  representing  the  product  of  all  the 
poles  except  the  one  for  which  the  partial  frac¬ 
tion  expansion  coefficient  is  being  sought  (and 
its  conjugate  if  complex).  The  odd  locations 
contain  the  coefficients  of  the  polynomial,  and 
the  even  locations  contain  the  corresponding 
power  of  s.  All  coefficients  are  real. 

2)  NP2  Number  of  occupied  locations  in  the  ?2 

array. 

8.  Subroutine  EVALU8  (P,  NP,  R,  V,  ZF) 
a.  Purpose 

This  routine  evaluates  a  polynomial  at  the  pole  for  which  the  par¬ 
tial  fraction  expansion  coefficient  is  being  sought. 
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b.  Input 


1)  P:  A  double  precision  array  containing  the  poly¬ 
nomial  coefficients  in  descending  order  with  the 
coefficient  (real  part,  imaginary  part)  and  power 
of  s  stored  in  separate  back-to-back  locations. 

2)  NP:  Number  of  occupied  locations  in  array  P. 

3)  R;  A  double  precision  array  containing  the  real 
part  of  the  pole  to  be  evaluated  in  the  first 
location  and  the  imaginary  part  in  the  next  loca¬ 
tion. 


c.  Output 

1)  V:  A  double  precision  array  containing  the  value 
of  the  polynomial  after  the  pole  of  Interest  is 
evaluated.  The  real  part  is  in  the  first  loca¬ 
tion  and  the  imaginary  part  in  the  second. 

2)  ZF:  A  scale  factor  that  is  used  in  the  routine 

to  maintain  numerical  accuracy  for  large  prod¬ 
ucts. 


9.  Subroutine  MULT IP  (Cl,  NT1 ,  C2,  NT 2 ,  C3 ,  NT3 ,  N) 


a.  Purpose 


This  routine  multiplies  two  polynomials  and  then  calls  the  subrou¬ 
tine  SIMPLE  to  combine  the  coefficients  with  like  powers  of  s. 


b.  Input 

1)  N:  An  integer  that  specifies  which  format  the 

polynomials  are  in.  A  two  corresponds  to  real, 
coefficients  with  two  locations  necessary  for 
each  polynomial  term.  A  three  corresponds  to 
real  and  imaginary  coefficients  with  three  loca¬ 
tions  necessary  for  each  polynomial  term. 

2)  Cl:  A  double  precision  array  containing  a  poly¬ 
nomial  in  the  format  specified  by  N. 

3)  NT1:  Number  of  occupied  elements  in  array  Cl. 
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4)  C2:  A  double  precision  array  containing  a  second 

polynomial  in  the  format  specified  by  N. 


5)  NT2:  Number  of  occupied  elements  in  array  C2. 
c.  Output 

1)  C3:  A  double  precision  array  containing  the 

product  of  the  polynomials  Cl  and  C2  in  the  format 
specified  by  N. 


2)  NT3:#)i)S:  Number  of  occupied  elements  in  array  C3. 


10.  Subroutine  GETPOL  (NR,  NRN,  A,  NA) 


a.  Purpose 


This  routine  takes  a  set  of  roots  and  multiplies  them  to  form  a 
polynomial. 

b.  Input 

1)  RN:  A  double  precision  array  containing  the 

roots  to  be  multiplied  together.  The  roots  are 
stored  real  part,  then  imaginary  part.  In  loca¬ 
tion  RN(2*NRN+1),  a  scale  factor  is  stored  that 
multiplies  the  polynomial. 

2)  NRN :  Number  of  roots  to  be  multiplied  together, 
b.  Output 

1)  A:  A  double  precision  array  with  coefficients  of 
the  polynomial  times  the  scale  factor  in  the  odd 
locations  in  descending  order  and  corresponding 
powers  of  s  located  in  the  even  locations. 

2)  NA:  Number  of  locations  used  in  array  A. 


1 1 .  Subroutine  ADD  (Cl,  NT1 ,  C2,  NT2,  C3,  NT3,  M) 


a.  Purpose 


The  routine  adds  polynomials  Cl  and  C2  together  and  places  the  sum 
in  array  C3.  The  polynomials  are  first  placed  sequentially  in  array  C3 
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and  then  the  subroutine  SIMPLE  is  called  to  add  coefficients  of  like 
powers  of  s.  The  dimension  of  array  C3  must  equal  the  combined  dimen¬ 
sions  of  Cl  and  C2. 

b.  Input 

1)  M:  An  integer  constant  specifying,  as  in  MULTIP, 
which  format  the  polynomial  coefficients  are  in. 

2)  Cl:  A  double  precision  array  containing  poly¬ 

nomial  coefficients  and  corresponding  powers  of  s 
in  the  format  specified  by  M. 

3)  NT1 :  Number  of  occupied  elements  in  array  Cl. 

4)  C2:  A  double  precision  array  containing  a  second 
polynomial  in  the  format  specified  by  M. 

5)  NT2:  Number  of  occupied  elements  in  array  C2. 

c.  Outputs 

1)  C3:  A  double  precision  array  containing  a  poly¬ 
nomial  which  is  the  sum  of  the  polynomials  Cl  and 
C2. 

2)  NT3:  Number  of  occupied  elements  in  array  C3. 

12.  Subroutine  SIMPLE  (P,  N,  K) 

a.  Purpose 

This  routine  combines  the  coefficients  of  a  polynomial  into  the 
least  number  of  coefficients  by  adding  together  the  coefficients  with 
like  powers  of  s. 
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b.  Input 


1)  K:  An  integer  constant  specifying,  as  in  MULTIP, 
the  format  of  the  polynomial  P, 

2)  P:  A  double  precision  array  containing  a  poly¬ 
nomial  in  the  format  specified  by  K. 

3)  N:  Number  of  occupied  elements  in  array  P. 

c.  Output 

1)  P:  Array  containing  least  number  of  coefficients 
necessary  to  specify  the  polynomial  read  in. 

2)  N:  Number  of  occupied  elements  in  the  output 

array  P. 

13 •  Subroutine  0RDER3  (P,  NP ,  (C) 

a.  Purpose 

This  routine  orders  the  polynomial  coefficients  into  descending 

powers  of  s. 

b.  Input 

1)  K:  An  integer  constant  specifying,  as  in  MULTIP, 
the  format  of  the  input  polynomial  P. 

2)  P:  A  double  precision  array  containing  the  input 
polynomial. 

3)  NP:  Number  of  occupied  locations  in  array  P. 

c.  Output 

1)  P:  Array  containing  the  polynomial  in  descending 
powers  of  s. 

2)  NP:  Number  of  occupied  locations  in  output 

array  P. 
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14.  Subroutine  CDEXP  (A,  B,  X,  Y) 


a.  Purpose 

This  routine  calculates  the  exponential  function  of  a  complex 
number  in  double  precision. 

b.  Input 

1)  A:  Real  part  (double  precision)  of  argument  of 

exponential  function. 

2)  B:  Imaginary  part  (double  precision)  of  argument 
of  exponential  function. 

c.  Output 

1)  X:  Real  part  (double  precision)  of  exponential 

function. 

2)  Y:  Imaginary  part  (double  precision)  of  expo¬ 

nential  function. 

15.  Subroutine  MULT  (A,  B,  C,  D,  X,  Y) 

a.  Purpose 

This  routine  multiplies  a  double  precision  complex  number  by  a 
double  precision  complex  number. 

b.  Input 

1)  A,  B:  Real  and  imaginary  parts  of  the  first 

number. 

2)  C,  D:  Real  and  imaginary  parts  of  second  number. 

c.  Output 

1)  X:  Real  part  of  the  complex  product. 

2)  Y:  Imaginary  part  of  the  complex  product. 
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16.  Subroutine  DIVI  (A,  B,  C,  D,  X,  Y) 


a.  Purpose 

This  routine  divides  a  double  precision  complex  number  by  a  double 
precision  complex  number. 


b.  Input 

1)  A,  B:  Real  and  imaginary  parts  of  the  first 

numbe  r . 

2)  C,  D:  Real  and  imaginary  parts  of  second  number. 


c.  Output 


17. 


18. 


1)  X:  Real  part  of  the  complex  division. 

2)  Y:  Imaginary  part  of  the  complex  division. 
Function  Fact  (n) 

This  is  a  function  subroutine  that  calculates  n! . 
Subroutine  POLYCO  (A,  B,  RR,  RI,  N) 


a.  Purpose 

The  purpose  of  this  routine  is  to  form  a  polynomial  from  a  set  of 
roots.  Both  real  and  imaginary  coefficients  are  calculated. 


b.  Input 

1)  RR:  Double  precision  array  containing  the  real 

part  of  each  root. 

2)  RI:  Double  precision  array  containing  the  imag¬ 
inary  part  of  each  root. 


3)  N:  Number  of  input  roots. 


c.  Output 

1)  A:  Double  precision  array  containing  the  real 

coefficients  of  the  polynomial. 

2)  B:  Double  precision  array  containing  the  Imagin¬ 
ary  part  of  the  polynomial  coefficients. 

19 •  Subroutine  ROOTS  (A,  B,  NN,  RR,  Rl) 

a.  Purpose 

This  subroutine  finds  the  roots  of  a  polynomial  with  complex  coef¬ 
ficients. 

b.  Input 

1)  A:  Double  precision  array  containing  the  real 

part  of  the  polynomial  coefficients  in  descending 
powers. 

2)  Double  precision  array  containing  the  imaginary 
part  of  the  polynomial  coefficients  in  descending 
powers . 

3)  NN:  Order  of  input  polynomial. 

c.  Output 

1)  RR:  Double  precision  array  containing  the  real 

part  of  each  root. 

2)  RI:  Double  precision  array  containing  the  imag¬ 
inary  part  of  each  root. 

F.  PROGRAM  VARIABLES  WITHIN  LABELED  COMMON 

Two  labeled  COMMON  blocks  are  used  in  the  main  program  ADVANZ,  and 
the  two  primary  subroutines  PARTFR  and  WPLN.  Variables  and  arrays  in 
the  main  program  and  the  two  primary  subroutines  share  the  same  storage 
locations  by  means  of  the  COMMON  statement.  These  variables  and  arrays 
are  stored  in  the  order  in  which  they  appear  in  the  block  specification. 
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The  COMMON  blocks  replace  the  subroutine  arguments  in  PARTFR  and  WPLN. 
This  arrangement  allows  DISCRET  to  be  easily  converted  to  an  overlay 
structure.  The  variables  and  arrays  located  within  these  COMMON  data 
blocks  are  outlined  below.  Each  variable  and  array  is  listed  indivi¬ 
dually  with  a  brief  description  of  its  purpose. 


COMMON/ ADVCZ/ RN , RD , NRTN , NRTD , MM , VR , VI , BR , BI , T , NH , AM , AO , TXFORM , CPLR 


Vari  able _ Purpose _  _ 

RN  A  double  precision  array  containing  the  zeros 

of  the  continuous  s-plane  transfer  function 
G(s).  The  zeros  are  stored  real  part  then 
imaginary  part  (required  storage  equals  two 
times  number  of  zeros). 

RD  A  double  precision  array  containing  the  poles 

of  G(s)  in  the  same  format  ae  the  zeros. 

NRTN  Number  of  zeros  in  G(s). 


NRTD 

MM 


VR,  VI 


BR,  BI 


T 


Number  of  poles  in  G(s). 

Integer  array  containing  multiplicity  of  poles 
corresponding  to  the  partial  fraction  expansion 
coefficients  located  in  arrays  VR  and  VI. 

Double  precision  arrays  containing  partial 
fraction  expansion  coefficients  (real  and  imag¬ 
inary)  for  the  poles  of  G(s)  including  those 
poles  Introduced  by  the  data  hold. 

Double  precision  arrays  containing  correspond¬ 
ing  poles  (reai  and  imaginary)  for  the  VR  and 
VI  arrays. 

Sampling  time  (sec). 


NH 


Integer  variables  used  in  the  data  hold  option. 


AO 


Gain  of  the  G(s)  transfer  function. 


TXFORM  Transform  option  variable  -  Z,  W,  or  WP. 

CPLR  Data  hold  option  variable  -  NON,  ZOH,  1ST,  2ND, 

or  SLE. 
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C0MM0N/T0TL12/CLNP0LY (51), CLDPOLY( 51 ) , CLZERQ( 50,2) 
CLP0LE( 50 , 2 ) , NCL2 , NCLP , CLK ,  CLNK , CLDK 


CLNPOLY 


CLDPOLY 


CLZERO 


CLPOLE 


NCLZ 


NCLP 


CLK 


Single  precision  array  containing  the  numerator 
polynomial  coefficients  for  the  discrete  trans¬ 
fer  function  G(z),  G(w),  or  G(w').  The  coeffi¬ 
cients  are  stored  sequentially  back-to-back, 
with  the  highest  order  coefficient  first  in  the 
array. 

Single  precision  array  containing  the  denomina¬ 
tor  polynomial  coefficients  for  the  discrete 
transfer  function  in  the  same  format  as  the 
numerator  array. 

Single  precision  array  containing  the  zeros  oi 
the  discrete  transfer  function.  The.  real  part 
of  the  nth  zero  is  stored  in  the  first  column 
(n,  1)  and  the  imaginary  part  in  the  second 
(n,2). 

Single  precision  array  containing  the  poles  of 
the  discrete  transfer  function  in  the  same 
format  as  the  zeros. 

Number  of  zeros  in  the  discrete  transfer  func¬ 
tion. 

Number  of  poles  in  the  discrete  transfer  func¬ 
tion. 

Total  gain  for  the  discrete  transfer  function 
(CLK  =  CLNK/ CLDK). 
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CLNK 


Numerator  gain  for  the  discrete  transfer  func¬ 
tion. 


CLDK  Denominator  gain  for  the  discrete  transfer 

function. 

G.  PROGRAM  OPERATING  INSTRUCTIONS 

The  example  in  Fig.  II  will  be  used  to  illustrate  the  input  and  out¬ 
put  data  structure  for  the  DISCRET  computer  program.  Input  data  items 
are  free-form  (f ree-format )  with  separators  rather  than  in  fixed-size 
fields.  The  two  exceptions  are  the  alphanumeric  inputs  which  select  the 
desired  data  hold  and  discrete  transform.  The  free-format  input  data 
consist  of  a  string  of  values  separated  by  one  or  more  blanks,  or  by  a 
comma  or  slash,  either  of  which  may  be  preceded  or  followed  by  any 
number  of  blanks.  A  line  boundary,  such  as  an  end  of  record  or  end  of 
card,  also  serves  as  a  value  separator  (Ref.  14). 

The  Input  is  divided  Into  three  main  blocks  of  data.  The  first 
block  contains  the  basic  parameters  that  define  the  s-plane  system  G(s) 
to  be  transformed.  These  data  are  placed  on  the  first  data  card  in  a 
free  format.  The  alphanumeric  code  for  the  data  hold  and  type  of  dis¬ 
crete  transform  are  inserted  on  data  cards  two  and  three  In  an  A3  and  A2 
format,  respectively.  The  final  block  of  data  is  again  free  format. 


T  ,T 

-  8ATs  G(s)M(s) 

Rt  J 


Figure  II.  Sampled  Continuous  System 
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This  input  starts  on  card  four  and  consists  of  the  zeros  and  poles  of 
the  continuous  s-plane  transfer  function  G(s).  The  required  input  data 
are  outlined  in  Table  2. 

The  arithmetic  sign  (i.e.,  positive  or  negative)  of  the  time  incre¬ 
ment  AT  selects  the  advanced  (positive  AT)  or  delayed  (negative  AT)  dis¬ 
crete  transforme  For  example,  for  a  sampling  period  of  T  *  1,0  and  a 
time  increment  of  0.3,  the  advanced  discrete  transform  is  specified  as 
AT  =  0.3  and  the  delayed  discrete  transform  as  AT  =  -0.3.  For  the 
standard  discrete  transform  the  time  increment  is  AT  =  0.0. 


TABLE  2. 

INPUT  DATA  FOR  DISCRET 

VARIABLE 

PURPOSE 

Data  Card  One  -  Free-Format 

NRTN 

Number  of  zeros  in  G(s) 

M 

Number  of  poles  in  G(s) 

T 

Sampling  time  (sec) 

A0 

G(s)  gain 

AM,  (AT) 

Time  increment  option: 

AT,,  -AT,  or  zero  (sec) 

Data  Cards  Two  and  Three  -  A3,  A2  Format 

CPLR 

Data  hold  option:  NON, 

ZOH,  1ST,  2ND,  or  SLE 

TXFORM 

Transform  option:  Z,  W, 
or  WP 

Data  Cards 

Four  to  nth  -  Free-Format 

RN(2*NRTN) 

Zeros  of  G(s):  real  part 
then  imaginary  part 

RD(2*M) 

Poles  of  G(s):  real  part 
then  imaginary  part 
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The  data  hold  and  transform  options  are  input  in  a  coded  alphanu¬ 
meric  format.  These  options  with  their  respective  alphanumeric  codes 
are  given  In  Table  3.  The  alphanumeric  input  code  for  the  data  hold  is 
placed  in  Columns  1-3  on  data  card  two  (left  justified).  The  discrete 
transform  code  appears  on  card  three  in  Columns  1-2  (left  justified). 


TABLE  3.  ALPHANUMERIC  CODES  FOR  DISCRET 


Data  Hold  Option 

Input  Code 

None 

NON 

Zero  order  hold 

ZOH 

First  order  hold 

1ST 

Second  order  hold 

2ND 

SI ewer 

SLE 

Transform  Option 

In nut  Code 

z  transform 

Z 

w  transform 

W 

w'  transform 

WP 

The  order  of  the  G(s)  transfer  function  must  be  equal  to  or  less 
than  50th,  Pole  multiplicity  up  to  and  Including  three  is  permitted. 
There  is  no  restriction  on  the  number  of  sets  of  repeated  poles.  The 
zeros  and  poles  of  G(s)  are  input  sequentially  in  a  free-format  (start¬ 
ing  on  data  card  four)  on  as  many  data  cards  as  is  necessary.  The  real 
and  imaginary  parts  are  separated  with  a  valid  separator  (i.e.,  a  comma, 
a  slash,  or  one  or  more  blanks).  The  zeros  are  given  first  followed  by 
the  poles.  For  real  roots,  0.0  must  be  input  for  the  imaginary  part. 

The  output  data  from  DISCRET  can  be  divided  into  two  main  sections. 
The  first  section  deals  with  the  input  parameters  for  the  s-plane 
continuous  system  C(s)  and  those  parameters  that  define  the  discrete 


transformation.  The  second  section  contains  che  results  of  the  trans¬ 
formation  process  -  the  discrete  transfer  function  G(z),  G(w),  or 

G(w').  The  program  first  prints  the  transform  options  that  have  been 
selected.  This  includes  the,  sampling  time,  data  hold,  type  of  discrete 
transform,  and  the  time  increment.  This  is  followed  with  a  list  of  the 
s-plane  zeros  and  poles  for  G(s)  including  those  introduced  by  the  data 
hold.  The  partial,  fraction  expansion  coefficients  for  G(s)  and  its  data 
hold  are  printed  next.  The  program  then  outputs  the  numerator  and 
denominator  polynomials  and  the  zeros  and  poles  for  the  discrete 
transfer  function. 

Table  4  contains  nine  sets  of  input  data  in  card  image  format.  The 
sampled  s-plane  systems  for  these  examples  are  depicted  in  Fig.  12.  The 
discrete  transfer  functions  for  each  of  these  systems  are  given  in 
Eqs.  73-75. 
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These  transfer  functions  are  calculated  by  the  DISCRET  computer  program. 

The  output  for  each  data  set  in  Table  4  Is  shown  in  Figs.  13-21. 
The  first  three  sets  of  data  calculate  the  standard  z-,  w-,  and  w'-  dis¬ 
crete  transforms  using  a  zero  order  hold  and  a  sampling  period  of 
T  ■=  0.1  (i.e.,  TXFORM  =  Z,  W,  and  WP;  CPLR  *  ZOH;  and  AT  «  0.0).  The 
output  for  these  three  examples  is  shown  in  Figs.  13-15.  In  data 
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sets  4-6,  the  advanced  z--,  w-,  and  w'-transforms  are  calculated  for  the 
second  system  given  by  Eq.  74.  A  ZQH  is  used  with  a  sampling  period  of 
T  =  0.04  and  a  time  advance  of  AT  “  0.004  seconds.  These  advanced  dis¬ 
crete  transfer  functions  appear  in  Figs.  16-18.  The  last  three  sets  of 
data  use  the  same  s-plane  system  G(s)  with  a  slewer  data  hold  and  a  time 
delay  of  AT  =  -0.004  (Eq.  75).  Figures  19-21  contain  the  delayed  dis¬ 
crete  transforms  for  these  inputs. 
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Figure  16a.  DISCRET  Output  for  Data  Set 
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Figure  16b.  DISCRET  Output  for  Data  Set 
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igure  17a.  DISCRET  Output  for  Data  Set 


Figure  18a.  DISCRET  Output  for  Data  Set 


caEfftciorrs  or  nunddatm  in 


Figure  18b.  DISCRET  Output  for  Data  Set 
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Figure  19a.  DISCRET  Output  for  Data  Set 


19b.  DISCRET  Output  for  Data  Set 


COEFFICIENTS  OF  NMBftTOft  IN 


Figure  21a.  DISCRET  Output  for  Data  Set 
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SECTION  V 


DESCRIPTION  OF  TXCONV  COMPUTER  PROGRAM 


A.  INTRODUCTION 

The  TXCONV  computer  program  calculates  a  low-rate  discrete  transform 
from  a  given  high-rate  discrete  transform.  The  input  to  TXCONV  is  a 

high-rate  transfer  function  in  the  z~,  w-,  or  w'-plane  and  the  output  a 
low-rate  transfer  function  in  the  z-,  w-,  or  w'-plane.  The  general 
transform  conversion  is  given  by: 

CT(z)  -  [cT/ffl(zp)]T  (76) 

CT(w)  =  [cT/m(Wp)]T  (77) 

CT(v')  =  [CT/m(Wp)]T  (78) 

The  superscript  designates  the  sampling  interval  used  to  form  the  dis- 

•T»  /  rn 

crete  transform.  For  example,  the  hlgh-rate  zp-transf orm  C  /m(rp)  is 
first  calculated  with  respect  to  a  T/m  sampling  period.  The  low-rate 
z-transform  of  this  high-rate  transform  is  then  taken  with  respect  to  a 
sampling  interval  of  T  seconds.  The  result  is  a  low- rate  z-transform 
CT(z) . 

The  transforms  in  Eqs.  76-68  are  generated  when  the  output  of  a 
system  is  sampled  at  a  lower  rate  than  the  input  (see  Section  III,  Sub¬ 
section  D).  An  open-loop  example  of  this  situation  is  depicted  in 
Fig.  22.  The  output  from  the  physical  sampler  T  in  Fig.  22  is  expressed 
as 


[GT/m]T  =  [GT/mRT/m]T 


(79) 
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R  RT/m 

c  /  c™/  [cT/m] T 

G 

T/m 

T/m  T 

( phantom) 

Figure  22.  High-Rate  Input /Low-Rate  Output  Sampling 

The  T/m  phantom  (fictitious)  sampler  in  the  output  facilitates  the  for¬ 
mation  of  the  GT/mgT/m  pro<juct  using  a  common  definition  of  the  trans¬ 
form  variable  (e.g.,  z^  ~  es^m).  This  mathematical  convenience  is 
valid  since  the  actual  output  sampler  T  simply  rejects  all  unwanted 
samples  from  the  phantom  sampler  T/m.  To  evaluate  Eq.  79,  the  procedure 
is  to  obtain  the  high-rate  T/m  transform  of  RT/,m  and  multiply  it  by  the 
high-rate  T/m  transform  of  GT//m.  This  high-rate  discrete  transfer  func¬ 
tion  product  is  the  required  input  to  the  TXCONV  computer  program.  The 
output  from  TXCONV  is  a  low-rate  discrete  transfer  function  defined  for 
a  T  sampling  period. 

TXCONV  is  written  in  FORTRAN  for  the  Control  Data  Corporation  (CDC) 
CYBER  175  series  computer.  The  program  can  handle  pole  multiplicity  up 
to  three  and  system  order  up  to  50th  (the  system  order  is  variable  and 
can  be  easily  changed).  There  is  no  restriction  on  the  number  of  sets 
of  repeated  poles.  Double  precision  arithmetic  is  used  throughout  the 
program. 

B.  PROGRAM  OPTIONS 

TXCONV  implements  the  conversion  of  a  h.igh-rate  discrete  transform 
to  a  low-rate  discrete  transfoim.  The  program  accepts  the  zeros  and 
poles  of  a  high-rate  discrete  transfer  function  in  the  z-,  w-,  or 

w'-plane  and  outputs  a  corresponding  low-rate  discrete  transfer  func¬ 
tion.  The  five  available  input  options  are  described  below. 


■  I 

:  i 


1 
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1 .  Z  Option 


The  program  input  consists  of  the  high-rate  zeros  and  poles  in  the 
z-plane  and  the  low-rate  (T)  and  high-rate  (T/ra)  sampling  periods. 
Prior  to  executing  the  transform  conversion,  the  high-rate  z-plane 
numerator  and  denominator  polynomials  are  transformed  to  the  w'-plane 
using  the  bilinear  transformation  zp  =  [ (2m/T)  4-  wp]/[(2m/T)  -  wp].  The 
w'-plane  denominator  is  then  rooted  to  obtain  the  poles  used  in  the 
residue  calculations.  In  this  option,  all  calculations  are  carried  out 
in  the  w'-plane  to  minimize  the  numerical  round-off  errors.  This  is 
necessary  since  the  poles  of  a  z-plane  function  tend  to  migrate  towards 
the  unit  circle  (i.e.,  z  **  1)  as  the  sample  rate  is  increased.  This  can 
introduce  numerical  errors  in  the  residue  computation.  These  inherent 
errors  can  be  minimized  by  performing  all  possible  calculations  in  the 
w'-plane  where  the  poles  are  more  reasonably  separated.  The  resulting 
low-rate  w'-plane  transfer  function  is  then  transformed  back  into  the 
z-plane  using  the  bilinear  transformation  w'  =  (2/T)(z  -  l)/(z  +  1). 
The  output  for  this  option  is  a  low-tale  discrete  transfer  function  in 
the  z-plane. 

2.  W  Option 

The  program  input  is  in  the  w-plane.  All  numerical  calculations  are 
carried  out  in  the  w-plane.  The  output  is  a  low-rate  discrete  transfer 
function  in  the  w-plane. 

3.  WF  Option 


The  program  input  is  in  the  w'-plane.  All  numerical  calculations 
are  carried  out  in  the  w'-plane.  The  output  is  a  low-rate  discrete 
transfer  function  in  the  w'-plane. 

4.  ZR  Option 

This  option  is  the  same  as  the  Z  option  except  that  the  w'-plane 
high-rate  poles  used  In  the  residue  calculations  are  obtained  by  direct 
transformation  of  the  input  z-plane  poles.  That  is,  the  high-rate 
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w'-plane  denominator  is  not  rooted  to  obtain  these  poles  as  is  done  in 
the  2  option.  The  ZR  option  avoids  the  numerical  errors  that  may  occur 
when  rooting  a  polynomial. 

5.  ZT  Option 

The  program  input  is  in  the  z-plane.  All  numerical  calculations  are 
carried  out  in  the  z-plane.  The  output  is  a  low-rate  discrete  transfer 
function  in  the  z-plane.  This  option  is  limited  to  simple  poles  (i.e., 
pole  multiplicity  equal  to  one). 

C.  TRANSFORMATION  EXPRESSIONS 

The  transformation  expressions  mechanized  in  the  TXCONV  computer 
program  are  given  below.  These  expressions  transform  a  high-rate  dis¬ 
crete  transfer  function  to  a  low-rate  discrete  transfer  function.  The 
mathematical  derivation  for  these  discrete  transformations  is  presented 
in  Section  III. 

1.  w  and  w'  Plane  Transformations 


residues 


CT/m(wp)  j 
An  +  v„  1  -  X 


A_  -  w_ 

w  “Poles  of  CT/m(e'  )/(An+w_) 


AP  +  wp  A  +  i/l 

Ac  -  “p  A  -  w  j 


Alternate  expressions  are  given  by 


CT(w)  «  J]  residues 


2Ap(A  +  w)[N(wp)/D*(wp)][l/(l  +  Yra 


v  +  A[ (1  -  Ym)/(1  +  Yra)l 


Vp-Foles  of  CT'/m(wp)/(Ap+Wp) 


Cr(w)  -  £  residues  2a  (A  +  v.)  j  — - 6 - —  j  (8. 

k  ^  I  D  (w  ) [ (w  +  A)  +  (w  -  A)Ym]  ,  r/m 

■  '  f  ’  Wp“Poles  of  C1'  <Wp)/(Ap-Hy 
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where 


C^’(w)  ■  N(w)/D(w) 


CT/m(wp)  -  N(Wp )/D(wp ) 


D* (wp )  -  D(wp)[(Ap  +  wp )(Ap  -  Wp ) j 


V  ■  (Ap  +  wp)/ (Ap  -  wp)  (87) 

These  transform  expressions  are  applicable  to  either  the  w-  or  w'-plane. 

For  the  w-plane,  the  transform  variables  are  w  and  wp  with  A  =  Ap  *  1. 

In  the  w'-piane,  the  transform  variables  become  w'  and  wp  with  A  -  2/T 

and  Ap  =  2m/T.  In  Eqs.  80-87  the  following  definitions  apply: 

m  =  Ratio  of  high-to-low  sampling,  (T)/(T/m) 


w,  w'  =  Low-rate  transform  variables 


Wp,  Wp  =  High-rate  transform  variables 


z  -  1 
z  +  1 


Zp  -  1 

#P  = 


2  z  -  1 
T  z  +  1 


2  zp  -  1 

(T/m)  zp  +  1 


,  zp  =  e«T/m 


2.  z-Plane  Transformations 


T/  x  V'  CT/®(Z  ) 

CMz)  =  resioues - - — 

k  zp  z  -  z?J 


Zp=Poles  of  C^/ra(Zp )/zp 


Alternate  expressions  are  given  by: 


T/  x  V*  I  N(z  )  ) 

Cx(z)  “  2~*  residues  z  r— — - - — ? - > 

k  | zD*(zp)  -  ZpO*(zp) j 


Zp»Poles  of  CT/m(Zp)/Zp 


CT(z) 


where 


-  E  residues  z  < 

k  ( 


N(zp)/D*(Zp)| 
2  ~  ZP  1 


Zp=Pcles  of  C^/m(zp)/zp 


CT(z)  =  N(z)/D(z) 


CT/mt^p)  N(zp)/D(zp) 


D  (zp)  —  ZpD(zp) 


=  Ratio  of  high-to-low  samplings  (T)/(T/m) 


sT 

z  =  Low-rate  transform  variable  (z  =  e  ) 


Zp  =  High-rate  transform  variable  (zp  =  c 


sT/m) 


D.  GENERAL  RESIDUE  CALCULATION 


Consider  the  general  partial  fraction  expansion  of  a  z-plane  func- 
tion  F(z)  for  3,  polo  with  multiplicity  ecjuai  to  ri . 


n _  +  An-1 _ ) _ An-2  ^ 


P  /  z  )  =  it  =  _ _ _ _ _  _ _ _ _ ^  _ _ 

D  (z  4-  a)n  (z  +  a)n_l  (z  +  a)n-'- 


An-3 

+ - ...  +  *  •  •  +  Remaining  poles  of  F(z) 

(z  +  a)n"3 
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The  partial  fraction  expansion  coefficients  are  given  by  "he  following 
equations  evaluated  at  z  =  -a: 


(n)  ! N 


D 


(95) 


An-l  - 


[  (n+1 )  ! / 1  !  ] N C ' )  1  - 


0>+1)d(’/ 


(96) 


|(n+2)  </2<  ]N'‘,)  -  AnD<'>"  "  -  (n+2)A„_1D(  ’ 


sr.+l 


(n+1)  (n+2)l)C)' 


(97) 


[  (n+3)  1/3!  J N < '  >  5  -  AT1P<2ri+3  -  (n+3)  A^jDC  ')n+'4  -  (n+2)(n+3)An_2D( ') 


n+2 


An-3 


\  n+l 


(n+1) (n+2) (n+3)D 


on 


(98) 


In  Eqs.  95-98,  (')n  is  defined  as  the  nth  derivative  with  respect  to  the 
transform  variable  z.  These  equations  are  completely  general  and  pro¬ 
vide  the  partial  fraction  expansion  coefficients  for  a  pole  with  multi¬ 
plicity  up  to  and  including  n  =  4.  The  coefficient  for  a  pole  with 
multiplicity  equal  to  n  =  1  is  given  by: 


A) 


_ N _ 

d(')1 


z=-a 


(99) 


For  n  =  2,  the  partial  fraction  expansion  coefficients  can  be  obtained 
from 


Al 


bN(  >^)1d(  ')2  -  2ND(  ')3 
3d(')2D(')2 


z=-a 


(100) 


A2 


(101) 
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And  for  n  *  3 


120n(')2dO3dO3  -  UNpO-V')5  -  SON^'^pO'V')4  +  I5H('^D( 

40d(')3d(')3d(')3 


(102) 


A 2  “ 


24n(')1d(')3  -  6ND(')4 


4D(')3D(')3 


(103) 


c=-a 


A3 


6N 


D(')3!z=-a 


(104) 


In  Eqs.  95-104,  the  pole  that  fs  being  evaluated  (z  +  a)n  is  not 
explicitly  factored  out  of  the  denominator  polynomial  D  prior  to  taking 
the  required  derivatives  or  prior  to  the  evaluation  at  z  =  -a.  The 
derivation  of  these  equations  follows  the  procedure  presented  in  Ref.  1. 
This  method  employs  L'Hopital's  rule  to  eliminate  the  indeterminate 
forms  that  result  (see  Ref.  1,  Appendix  D  for  details).  The  general 
procedure  is  to  take,  consecutive  derivatives  of  the  numerator  and  denom¬ 
inator  polynomials  until  a  determinate  form  Is  obtained.  For  example, 
in  the  n  =  1  case,  we  have 

N  Ai 

F(z)  -  ~  =  (z~+~a7  +  *’*  +  Remaining  poles  of  F(z)  (105) 


Then 


A1 


(z  +  a)N 
D 


z=-a 


(106) 


If  Eq.  106  is  evaluated  a  z  =  -a  without  first  explicitly  factoring  out 
the  (z  +  a)  factor  in  the  denominator  and  cancelling  it  with  the  numera¬ 
tor  factor  (z  4-  a),  an  indeterminate  form  (0/0)  will  result.  Applying 
L'Hopital's  rule  (i.e.,  taking  separate  derivatives  of  the  numerator  and 
denominator  with  respect  co  z)  gives 
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(z  +  a)N '  +  N 
D7" 


007) 


Notice  that  for  any  pole  with  a  multiplicity  of  n  =  1,  the  evaluation  of 
the  first  derivative  of  the  denominator  produces  a  finite  result.  That 


D'l  „  t  0.0 

1  —  3 


(108) 


This  can  be  stated  in  more  general  terms  for  any  pole  with  multiplicity 
equal  to  n  as 


D(')k  __  ]f  0.0  ,  (z  +  a)n 


k  >  n 


(109) 


D(')k 


-  0.0 

lz=-a 


(z  +  a)n 


k  <  n 


(110) 


where,  again,  (';  is  defined  as  the  kth  derivative  of  the  denomina¬ 


tor  D. 


For  example,  if 


N_  -  (  +  2)(z  +  3) 

D  (z+5)(z+10) 


(z  +5)  (z  +  10) 


F(z ) 


(111) 


Then 

D  =  z2  +  15z  4-  50 
D'  =  2z  4-  15 


and 


Al 


N 

D' 


z=-10 


z^  +  5z  +  6 

2z  4-  15  z=— iO 


-13.2 


Bl 


N_ 

D' 


z=-5 


z^  +  5z  4-  6 
2z  -  15 


z=-5 


1.2 


(112) 

(113) 


(114) 


(115) 


The  evaluation  of  the  residues  in  the  TXCONV  computer  program  is 
accomplished  using  Eqs.  99,  100,  and  102.  It  is  recognized  that  the 

residue  for  a  pole  with  multiplicity  equal  to  n  is  given  by  the  partial 
fraction  expansion  coefficient  associated  wilii  the  (z  4-  a)  terra.  Apply¬ 
ing  Eqs.  99,  100,  and  102  to  Eqs.  80  and  88  results  in  closed-form  solu¬ 
tions  which  are  functions  of  the  separate  derivatives  of  the  numerator 
and  denominator  polynomials  for  the  given  high-rate  discrete  transfer 
function.  The  actual  expressions  mechanized  in  TXCONV  are  developed  in 
the  next  subsection. 

E.  MECHANIZATION  SCHEME 

The  transformation  expression  in  the  z-plane  is  only  mechanized  for 
poles  with  multiplicity  equal  to  one  (ZT  option).  This  transformation 
is  coded  in  subroutine  ZMULT1.  As  explained  previously,  for  pole  multi¬ 
plicity  greater  than  one  (n  >  1),  the  input  high-rate  z-plane  transfer 
function  is  first  transformed  to  the  w'-plaue  prior  to  performing  the 
numerical  calculations  (Z  and  ZR  options).  This  procedure  improves  the 
accuracy  of  the  results  by  reducing  to  a  minimum  the  errors  introduced 
by  numerical  round-off.  The  z-plaue  convesion  for  option  ZT  is  formed 
by  applying  Eq.  99  to  Eq.  88.  The  resul  is  Eq.  116  (see  Subsections  C 
and  D  for  definition  of  terras) . 
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CT(z) 


(H6) 


T.  residues 
k 


M(zp )/D*(zp) ' | 

z  “  zp  j 


Zp=Poles  of  D*(Zp) 


where 


zpD(zp) 

017) 

(T)/(T/m) 

(118) 

esT  ,  zp  =  esT/m 

(119) 

For  all  options  except  ZT  (i.e.,  Z,  W,  WP,  and  ZR),  separate  mechan¬ 
ization  expressions  are  used  for  poles  with  multiplicity  equal  to  one, 
two,  and  three.  These  expressions  are.  coded  in  the  subroutines  WMuLTi , 
WMULT2,  and  WMULT3  for  n  =  1,  2,  and  3,  respectively.  Equations  99, 
100,  and  102  are  individually  applied  to  Eq.  80  to  form  these  expres¬ 
sions.  The  resulting  transform  conversion  equations  (Eqs.  120,  121,  and 
124)  are  applicable  to  either  the  w-  or  w'-plane.  The  definition  of 
terms  for  w  or  w'  implementation  is  given  in  Subsections  C  and  D.  The 
specific  expressions  that  are  mechanized  are  outlined  below. 

Multiplicity  n  =  1 


C1(w) 


)  ,  residues  2A(1(A  +  W) 
k 


N(wp)/D*(wp)'[l/(1  +  Yra)J 
w  +  TOLE 


wp=roles  of  D(up) (Ap+Wp) 


(120) 


8  7 


Multiplicity  of  n  =  2 


where 


^residues  2A  (A  +  W)  -N— +- 
k  p  u  +  (POLE)2 


vp-Poletf  of  D(wp)  (Ap-HWp) 


(121) 


2N 1  (w  ) 


2N(v  )D*f"  (w  )  2raN(wJYm''1Y' 


D*"(wp)  (1  +  Yra)  3(D*"(«p))2(l  +  Y™)  D*"(wp)(i  +  Ym)  2 


(122) 


2N’  (w  )  (POLE)  2N(wn)D*"'  (w  )  (POLE)  2mAN(wn)Yn'  V 

N2  - E - 2 — — - E -  +  - E - 

D*"(wp) (X  +  Yra)  3(Urt"(wp))2(l  +  Ym)  D*"(wp) (1  +  Yra)2 


(123) 


Multiplicity  of  n  -  3 


wher  e 


) (  2Ap(A  f- 


{(HI  +  N2)[w2  +  2  (POLE)w  +  (POLE)2] 

+  (N3  +  N4)  [ w2  +  (POLE  -  A)w  -  A(POLE)] 

-f  (N5)  [w2  -  (2A)w  +  A2) _ __ _ 

V  u  +  (POLE)3 


(124") 


Wp«Poles  of 
D(wp)  (Ap+Wp) 


3N"(wn) 


X.5N'(w0)D*""(wd) 


D*'"  («_)(!  +  Ym)  (D*’"  (w  ))2(i  +  Ym) 
P  P 


(125) 


■  375N(wp)(D*""(wp))2  0.3N(wp)D*  '""(wp) 

(D*"'  (wp))3(l  +  Ym)  (D*m  (wp))2(l  +  Yra) 


(126) 


3mN (v  ) Ym'1Y ' D*" " 


'  (w  ) Ym~1Y  *  X , 5mH (w  ) Ym_1Y ' D*""(w  ) 


(D*’"  (wp))2(l  +  f™)2  D*"‘  (wp)  (1  +  Ym)2  (E*'"  (wp))2(X  +  Y»)2 


(127) 


3m(m  -  l)H(w  )Yn_2(Y')2  3mN(wp)Y,n~lY" 
D*"'  (wp)  (X  +  Ym) 2  D*"‘  (vp)(l  +  Ym)2 


(128) 


6m2N(wp)  (Y"1”1)  (Y')2 
»*'"  (wp)(l  +  Y“)3 


(129) 
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For  Eqs .  120-129,  the  following  relationships  and  definitions  apply: 


POLE  =>  A[  ( 1  -  Ym)/(1  +  Ym)j 

(130) 

Y  -  (Ap  +  wp)/(Ap  -  wp) 

<130 

Y  =  2Ap/ (Ap  —  Wp)^ 

(132) 

Y"  =•  4Ap/(Ap  -  wp)3 

(133) 

D*(wp /  =  D(wp)(Ap  +  wp ) ( Ap  -  wp) 

(134) 

m  =  (T)/(T/m) 

(135) 

F.  ROOT  CONVERSION  BETWEEN  s,  z,  w,  or  w'  PLANES 

The  conversion  of  the  poles  between  the  s-,  z-,  w-,  or  w'-plane  is 
mechanized  in  the  subroutine  SZWR00T.  The  algebraic  relationships 
implemented  Ir  SZWR00T  allow  direct  transformation  of  the  poles  from  one 
complex  plane  to  another.  Direct  conversion  of  the  zeros  on  a  one-for- 
one  basis  is  not  normally  possible.  One  exception  is  the  conversion 
between  the  w-  and  w'-plane  where  both  the  zeros  and  poles  transform 
directly  (w'  =  2w/T).  It  is  also  possible  to  directly  transform  the 
zeros  and  poles  from  the  w-  or  w'-plane  to  the  z-plane.  However,  the 

J  —  -  -  -  C-  L  —  ’  *_  - 11..  /  _  -  1  _ -  -  -  -  -  '  _  1  _  _ 

mveibe  uub  UHL  uuuudiiy  lluc  \  x  •  t:  •  ,  t  “picnic  t_u  w—  vi  w  -pxaui; 

direct  conversion  of  zeros). 

The  complete  conversion  between  any  of  the  complex  planes  noted 
above  can  be  implemented  by  transforming  the  complex  function  instead  of 
its  zeros  and  poles.  Although,  in  some  cases,  physically  unrealizable 
functions  may  result.  The  general  procedure  is  to  transform  the  numera¬ 
tor  and  denominator  polynomials  (cr  partial  fraction  expansion  terms) 
and  then  factor  the  results  to  obtain  the  transformed  zeros  and  poles. 
This  procedure  is  normally  less  desirable  than  direct  conversion  of  the 
zeros  and  poles,  since  the  particular  mathematical  calculations  (e.g., 
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numerical  factoring;  multiplication  and  addition  of  polynomials)  may 
introduce  additional  numerical  errors. 

The  subroutine  SZWROOT  is  used  in  the  TXCONV  computer  program  to 
transform  the  poles  between  the  z-  and  w '-planes.  The  other  conversions 
available  in  SZWROOT  are  not  utilized.  Conversion  of  the  poles  between 
the  z-  and  w '-planes  occur  only  in  the  ZR  option.  The  algebraic  conver¬ 
sion  relationships  implemented  in  SZWROOT  are  outlined  below.  The 
following  definitions  apply  to  the  root  conversion  equations: 


z  ^  esT  =  !_+_5i  *  (2/T)  +  w' 

1  -  w  (2/T)  ~  w' 


(136) 


»  ■  Hr  ■  <T/2)”' 


(137) 


w 


2  z  -  1 
T  2  +  1 


(2/T)w 


(138) 


s  -  X6  +  jYs  ,  z  »  Xz  +  jYz  (139) 


w  =  Xw  +  jYw  ,  w'  =  XW'  +  jYw'  (140) 


3  -Plane  to  z-Plane 


X 


z 


eXsT  cos  YST 


(141) 


Yz  «  eX®T  sin  YST 


(142) 
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w-Plane  to  z-Plane 


1  ~  x£  -  y£ 

(1  -  XJ2  +  Y2 


(143) 


*z 


2Yl 


(1  -  Xw)2  +  Y2 


(144) 


w'-Plane  to  z-Plane 


x  =  1  -  [ (T/2)XW/ ] 2  -  [ (T/2)Yw<]2 

[l  -  (T/2)XW']2  +  [(T/2)Yw']2 

Y  =  2[(T/2)YW'] _ 

lz  [l  -  CT/2)XW- ] 2  +  [ (T/2)Yw']2 

z-Plane  to  s--Pl.ane 


(145) 


(146) 


Xs  -  ^  In  [xl  +  Y2z)  (147) 

Ys  =  i  tan'l  (yz/xz)  ( 148 > 


w^Plane  to  s-Plane 


Xs 


1  (1  +  xw)2  +  *2 

2T  "  (1  -  Xw)2  +  Y2 


(149) 


1  _  2Yw 

7r  tail  1 - y— 

T  1  ~  X2  ■ 


Lw 


(150) 


91 


w'-Plane  to  s-Plane 


1  t  l_J_  (T/2)XW^]2  +  [  (T/2)YW']2 
3  *  2T  ln  [i  _  (T / 2 ) Xw ' ] 2  +  [ (T/2) Yw' ] 2 


1  „  -!  2[(T/2)YwJ _ 

S  T  1  -  [(T/2)XW']?  -  [(T/2)YW']2 


s-Plane  to  w-Plane 


XoT  -XoT 
e  a  -  e  ® 


XCT  .  -XCT 


+  e  s  +  2  coa  YoT 


2  stn  YST 


v  m  Y  m 

e  s 1  +  e”*8A  +  2  cos  YRT 


z-Plane  to  w-Plane 


Xw 


xj  +  Y2  -  i 
(X2  +  l)2  +  Y2 


*w 


2YZ  _ 

(Xz  +  l)2  +  Y2 


(151) 


(152) 


(153) 


(154) 


(155) 


(156) 
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w'-Plano  Co  w-Plane 


(157) 


Y 


w 


T 

7  Yw' 


(158) 


s-Plane  Co  w'-Plane 


XST  _  ~XST 


T  XoT  ,  -XST  ,  •>  v  t 
e  s  +  e  s  -f-  2  cos  YST 
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2  sin  YST 


T  eXsT  +  e  XsT  +  2  cos  YST 


'  1  /Lf\\ 
i  x  w  J 
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G.  PROGRAM  STRUCTURE 


The  basic  structure  of  the  TXCONV  computer  program  contains  a  main 
program  TXCONV  and  seven  primary  subroutines  ZMULT1,  WMULT1,  WMULT2, 
WMULT3,  RES1,  RES2,  and  RES3.  The  program  modules  along  with  the 
20  supporting  subroutines  are  depicted  in  Fig.  23.  The  main  program 
TXCONV  first  calls  ZMULTl,  WMULT1,  WMULT2 ,  and  WMULT3  to  evaluate  the 
individual  residues  of  the  high-rate  poles.  Subroutines  RES1,  RES2,  and 
RES3  are  then  called  to  combine  these  residues  to  form  the  low-rate  dis¬ 
crete  transfer  function.  ZMULTL  calculates  the  residues  for  poles  with 
multiplicity  equal  to  one  (n  =  1).  This  subroutine  mechanizes  Eq.  116 
in  the  z-plane.  WMULT1,  WMULT2,  and  WMULT3  calculate  the  residue.s  for 
poles  with  multiplicity  equal  to  one,  two,  and  three  (n  =  1,  2,  3),  res¬ 
pectively.  These  subroutines  mechanize  Eqs.  120,  121,  and  124  in  the  w- 
or  w'-plane. 

The  TXCONV  computer  program  can  be  operated  in  a  standard  program- 
subroucitie  structure  (i.e.,  a  main  program  followed  by  its  subroutine) 
or  in  an  overlay  structure.  The  source  listing  for  TXCONV  in  Volume.  Ill 
is  set  up  in  a  program-subroutine  structure.  However,  this  listing  also 
includes  (in  the  comment  code)  the  required  changes  to  run  the  program 
in  an  overlay  structure.  Dividing  the  program  into  overlays  reduces  the 
amount  of  computer  memory  required  to  execute  the  program.  The  overlay 
code  is  highlighted  with  a  star  (*)  character  in  column  one  (which  makes 
the  code  inactive  comment).  Removing  this  code  from  comment  will  allow 
overlay  operation.  To  complete  this  turnover,  the  TXCONV  program  card 
and  subroutine  cards  for  ZMULTl,  WMULT1,  WMULT2,  WMULT3 ,  RES l,  RES2,  and 
RES3  must  be  deleted.  In  addition,  the  call  statements  for  these  sub¬ 
routines  (which  are  now  overlays)  must  also  be  deleted.  These  call 
statements  are  located  in  the  "Master  Do  Loop  for  Residues"  section  in 
the  main  program.  The  overlay  code  creates  a  main  overlay,  a  primary 
overlay,  and  seven  secondary  overlays.  This  overlay  structure  is  out¬ 
lined  below: 
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Figure  23.  TXCONV  Frogram  Structure 


Overlay 

Program  Name 

(TXCONV, 0,0,0) 

TXCONV 

(27,0) 

MAIN 

(27,1) 

ZMULTl 

(27,2) 

WMULT1 

(27,3) 

WMULT2 

(27, A) 

WMULT3 

(27,5) 

RESl 

(27,6) 

RES2 

(27,7) 

RES3 

Parameters  are  passed  between  the  main  program  and  the  seven  primary 
subroutines  entirely  in  a  common  data  structure.  This  lack  of  formal 
parameter  passing  via  subroutine  arguments  was  initiated  to  allow  the 
program  to  be  easily  converted  to  an  overlay  structure  in  the  TOTAL 
(Ref.  13)  computer  program  at  AFWAL/FIGC.  This  permits  interactive 
operation  of  TXCONV  as  an  option  in  TOTAL. 

H.  DESCRIPTION  OF  SUBROUTINES 

This  subsection  contains  a  brief  definition  of  the  routines  used  in 
the  TXCONV  computer  program.  A  detailed  description  of  each  routine  is 
documented  in  the  COMMENT  section  of  the  source  code.  This  comment  code 
defines  all  subroutine  arguments  and  the  major  internal  variables  in 
each  routine.  In  addition,  descriptive  comments  are  included  throughout 
each  routine.  The  present  source  code  (Volume  III)  is  set  up  to  handle 
50th  order  systems.  The  comment  code  in  the  main  program  TXCONV  lists 
the  required  changes  to  the  program  DIMENSION  statements  to  alter  the 
maximum  system  order  allowed.  The  purpose  of  each  variable  and  array  in 
labeled  COMMON  is  also  included  in  the  source  code  for  the  main  program. 

Program  TXCONV.  This  is  the  main  program  for  the  TXCONV  computer 
ptogram.  It  reads  the  input  from  data  cards  and  transfers  the  data  to 
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internal  variables  and  arrays  that  are  located  in  the  COMMON  data 
structure.  It  adds  appropriate  denominator  poles  (z  -  0,  w  -  -1,  or 
w'  =  2m/T)  to  the  residue  expressions.  The  call  statements  for  the 
seven  primary  subroutines  (ZMULTI,  WMULT1,  WMULT2,  WMULT3,  RES1,  RES2, 
and  RES3)  are  located  in  the  main  program.  The  final  formation  of  the 
output  low-rate  discrete  transfer  function  is  accomplished  in  this  pro¬ 
gram  module.  The  individual  low-rate  discrete  transfer  functions  from 
subroutines  RES1,  RES2,  and  RES3  are  combined  to  form  the  final  output. 

Subroutine  ZMULTI.  This  subroutine  calculates  the  residues  for 

poles  with  multiplicity  equal  to  one  (n  =  1).  ZMULTI  mechanizes  the 

z-plane  residue  expression  in  Eq.  116. 

Subroutine  WMULT1.  This  subroutine  calculates  the  residues  for 

poles  with  multiplicity  equal  to  one  (n  =  1).  WMULT1  mechanizes  the 

w-plane/w'-plane  residue  expression  in  Eq.  120. 

Subroutine  WMULT2.  This  subroutine  calculates  the  residues  for 

poles  with  multiplicity  equal  to  two  (n  =  2).  WMULT2  mechanizes  the 

w-plane/w'-plane  residue  expression  in  Eq.  121. 

Subroutine  WMULT3.  This  subroutine  calculates  the  residues  for 

poles  with  multiplicity  equal  to  three  (n  =  3).  WMULT3  mechanizes  the 

w-plane/w'-plane  residue  expression  in  Eq.  124. 

Subroutine  PJSl.  This  subroutine  forms  the  overall  low-rate  dis¬ 
crete  transfer  function  for  1st  order  poies.  The  n  =  1  (multiplicity 

equal  one)  residues  that  are  calculated  in  subroutines  ZMULTI  or  WMULT1 

are  combined  to  form  this  transfer  function. 

Subroutine  RES2.  This  subroutine  forms  the  overall  low-rate  dis¬ 
crete  transfer  function  for  2nd  order  poles.  The  n  -  2  (multiplicity 

equal  two)  residues  that  are  calculated  in  subroutine  WMULT2  are  com¬ 
bined  to  form  this  transfer  function. 

Subroutine  RES3.  This  subroutine  forms  the  overall  low-rate  dis¬ 
crete  transfer  function  for  3rd  order  poles.  The  n  =  3  (multiplicity 

equal  three)  residues  that  are  calculated  in  subroutine  WMULT3  are  com¬ 
bined  to  form  this  transfer  function. 
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Subroutine  COMPOLY.  Thin  subroutine  forms  a  polynomial  from  a  set 
of  roots.  Both  real  and  imaginary  coefficients  are  calculated. 

Subroutine  BILIN.  This  subroutine  performs  the  general  bilinear 
transformation  from  one  complex  plane  to  another. 

Subroutine  TERM.  This  subroutine  calculates  the  individual  terms 

used  in  the  bilinear  transformation  mechanized  in  subroutine  BILIN. 

Subroutine  WZBILIN.  This  subroutine  initiates  the  specific  bilinear 
transformation  from  the  w'-plane  to  the  z-plane.  The  actual  transforma¬ 
tion  is  carried  out  in  subroutine  BILIN. 

Subroutine  MULTIP.  This  subroutine  multiplies  two  polynomials. 

Subroutine  ADD.  This  subroutine  adds  two  polynomials. 

Subroutine  SIMPLE.  This  subroutine  simplifies  a  polynomial  by  add¬ 
ing  coefficients  of  like  powers. 

Subroutine  COEFF.  This  subroutine  adds  the  missing  power  terms  in  a 
polynomial  by  inserting  a  zero  coefficient  with  the  appropriate  power 
and  moving  the  original  terms  to  make  room  for  the  missing  terms. 

Subroutine  ORDER3 .  This  subroutine  orders  a  polynomial  in  descend¬ 
ing  powers. 

Subroutine  HOOTS.  This  subroutine  finds  the  roots  of  a  polynomial 
with  complex  coefficients. 

Subroutine  SZ-WROOT.  This  subroutine  performs  the  root  conversion 
between  the  s-,  z-,  w-,  and  w'-complex  planes. 

Subroutine  ORDPOLE.  This  subroutine  checks  for  multiple  poles  and 
stores  the  multiplicity  in  an  array.  The  extra  multiple  poles  are 
deleted  and  only  a  single  copy  of  each  pole  is  stored  in  the  output 
array. 

Subroutine  POLE.  This  subroutine  calculates  a  low-rate  pole  in  the 
z-plane  from  a  given  htgh-rate  pole  in  the  z-plane. 

Subroutine  WPOLE.  This  subroutine  calculates  a  low-rate  pole  in  the 
w-  or  w'-plane  from  a  given  high-rate  pole  in  the  w-  or  w'-plane. 
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Subroutine  DERIV3.  This  subroutine  takes  the  derivative  of  a  poly¬ 
nomial. 


Subroutine  EVALU3.  This  subroutine  evaluates  a  polynomial  for  a 
given  root. 

Subroutine  DIVI.  This  subroutines  divides  two  complex  numbers. 

Subroutine  MULT.  This  subroutine  multiplies  two  complex  numbers. 

Subroutine  DOLOOP .  This  subroutine  implements  a  standard  DO  LOOP  to 
transfer  one  array  Into  another  array. 

Subroutine  CANROQT.  This  subroutine  cancels  equal  zeros  and  poles 
according  to  a  specified  tolerance.  Separate  tolerances  are  provided 
for  the  real  and  imaginary  parts  of  each  root. 

I.  PROGRAM  OPERATING  INSTRUCTIONS 

This  subsection  presents  the  input  and  output  data  structure  for  the 
TXCONV  computer  program.  Most  input  data  are  in  free  format  with  sepa¬ 
rators  rather  than  in  fixed-size  fields.  The  one  exception  is  the 
alphanumeric  input  which  selects  the  transformation  option  (Z,  W,  WP, 
ZR,  or  ZT ) .  The  free-format  input  data  consist  of  a  string  of  values 
separated  by  one  or  more  blanks,  or  by  a  comma  or  slash,  either  of  which 
may  be  preceded  or  followed  by  any  number  of  blanks.  A  line  boundary, 
such  as  an  end  of  record  or  end  of  card,  also  serves  as  a  value 
separator  (Ref.  14). 

The  first  section  of  input  data  contains  the  parameters  that  define 
the  high-rate  discrete  transfer  function,  the  high-rate  sampling  period 
(TIN),  and  the  low-rate  sampling  period  (TOUT).  These  data  are  placed 
on  the  first  two  data  cards  in  a  free  format.  The  alphanumeric  code  for 
the  transformation  option  Is  placed  on  data  card  three  in  an  A2  format. 
(See  Subsection  B  for  an  explanation  of  _he  transformation  options.) 
The  remaining  data  are  again  free  format  and  consist  of  the  zeros  and 
poles  of  the  high-rate  transfer  function.  The  order  of  this  transfer 
function  (i.e.,  order  of  the  denominator  polynomial)  can  be  equal  to  or 
less  than  50  with  pole  multiplicity  up  to  and  including  three.  .Mu 
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zeros  and  poles  are  input  sequentially  in  a  free  format  (starting  on 

data  card  four)  on  as  many  data  cards  as  necessary.  The  real  and  imag¬ 
inary  p'-.rts  are  separated  with  a  valid  separator  (i.e.,  a  comma,  a 

slash,  or  one  or  more  blanks).  The  zeros  are  inserted  first  followed  by 

the  poles.  For  real  roots,  0.0  must  be  input  for  the  imaginary  part. 

The  required  input  data  are  outlined  below: 


Data  Card 


Data  Items 


1 

2 

3 

4 

nth 


GAIN,  NZEROS,  NPOLES 
TIN,  TOUT 

Transformation  option 
Zeros(real,imag) 
Poies(rea’  nag) 


Format 
Free  format 
Free  format 
A2 

Free  format 
Free  format 


The  following  definitions  apply  to  the  data  items  listed  above: 


GAIN  -  High-rate  transfer  function  gain 

NZEROS  -  Number  of  zeros 

NPOLES  -  Number  of  poles 

TIN  -  High-rate  sampling  period  (sec) 

TOUT  -  Low-rate  sampling  period  (sec) 

Transformation  option  -  Z,  W,  WP,  ZR,  or  ZT 

Zeros (real, iraag)  -  Real  and  imaginary  parts  of 
high-rate  zeros 

Poles ( real, imag)  -  Real  and  imaginary  parts  of 
high-rate  poles 


The  output  data  from  TXCONV  are  divided  into  two  main  sections.  The 
first  section  contains  a  listing  of  the  input  variables  associated  with 
the  high-rate  discrete  transfer  function.  These  variables  include  the 
ratio  of  the  sampling  periods  (RATIO);  the  high-rate  (TIN)  and  low-rate 
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id* 


(TOUT)  sampling  periods;  the  high-rate  transfer  function  gain  (GAIN); 
and  the  number  of  input  zeros  (NZEROS)  and  poles  (NPOI.ES).  The  program 
also  prints  the  numerator  and  denominator  polynomials  and  roots  for  the 
high-rate  discrete  transfer  function.  A  listing  of  the  high-rate  poles 
with  their  multiplicity  is  then  printed.  An  additional  high-rate  pole 
(required  by  the  transformation  process)  at  z  *  0.0,  w  =  1.0,  or 

w'  =  2/TIN  is  included  in  this  listing  (see  Section  III).  The  second 
section  of  output  data  deals  with  the  low-rate  discrete  transfer  func¬ 
tion.  The  numerator  and  denominator  polynomials  and  their  roots  are 
printed.  The  low-rate  transfer  function  gain  and  a  second  listing  of 
the  high-rate  and  low-rate  sampling  periods  are  also  given. 

The  fast-input/slow-output  sampled  system  in  Fig.  24  will  be  used  to 
illustrate  the  input  and  output  data  structure  for  TXCONV .  The  proce¬ 
dure  for  this  example  is  typical  for  closed-loop  systems  employing  fast- 
input/slow-output  sampling. 


Figure  24.  Fast-Input/Slow-Output 
Sampled  System 
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The  output  equation  for  Fig.  24  is  given  by 


C  -  (GiM2]eT/2  (165) 

where 

eT/2  .  rT/2  _  (G2M)T/2(G1M2E'r/'2)T  (166) 

To  solve  for  the  ET^2  signal  in  Eq.  166,  premultiply  by  G^M2  and  take 
the  T  transform  of  both  sides  of  the  resulting  equation  (or  sample  both 
sides  of  the  equation  at  a  T  interval).  The  result  is 

(g^E11/2)1  -  (g1M2RT/2)T  -  [(G1M2)(G2M)T/2]T(GiM2ET/2)T  (167) 

Rearranging  Eq.  167, 

(g1M2Et/2)T  -  j  1  +  [(G1M2j(G2M)T/2]Tj  1  (GiM2Rt/2)T  (168) 

and  substituting  Eq.  168  into  Eq.  166  produces  Eq.  169. 

ET/2  -  rT/2  _  (g2m)t/2  | I  +  [(G1M2)(G2M)T/2]J |  *  (G1M2Rt/2)T  (169) 

Finally,  substituting  Eq.  169  into  Eq.  165  gives  the  output  equations 
for  Fig.  24. 

C  «  (g1M2)RT/2  -  (G1M2)t/2(G2M)t/2  |i  +  [[G1M2)(G2M),r/2]Tj  1  (G1M2RT/2)T 

(170) 
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To  illustrate  the  operation  of  the  TXCONV  computer  program,  we  cal¬ 
culate  the  term 

[(GiM2)(G2M)T/2]T  (171) 

I,et  z  =>  es^2>  and  introduce  a  phantom  T/2  sampler  to  Eq.  171.  This 
mathematical  operation  is  depicted  in  Fig.  25.  This  step  is  valid  since 
the  T  output  sampler  simply  rejects  all  the  unwanted  samples  from  the 
T/2  sampler.  Equation  171  then  becomes 


[(GiM2)T/2(G2M)T/2]T 


(172) 


R  =  l 

mg2 

/ 

M2G, 

/  /  ^ 

T/2 

T/2  T 

(phantom) 


Figure  25.  Phantom  Sampler  Concept 


We  next  choose  a  low-rate  sampling  period  of  T  *  -ln(0,81).  This  gives 
a  set  of  convenient  numbers  when  Eq.  172  is  evaluated.  Inserting  the 
transforms  depicted  in  Fig.  24  into  Eq.  172  produces 


(GiM2)T/2 


e-1'/2 


(G2M)t/2 


e"T/2 


z  +  1 
z 


(173) 


z 


z 


and 


[(GlM2]T/2(G2M)T/2]T 


( 1  -  e~T/2)2(z  +  i) 
z(z  -  e“^/2j2 


,T 


2  qST/2 


0.01(2  +  1) 

z(z  -  0.9)2 


074) 


Equation  174  can  be  solved  by  calculating  the  residues  of  the  following 
expression  (see  Eqs.  36  and  88). 


1 

2irj 


(l  -  e“T/2)2(2jp^fJO  z  dz 

zp[zp  ”  e~^/2j2  (z  ~  zp)  *-p 


(175) 


The  residues  for  the  double  poles  at  zp 


0.0  and  zp 


e-T/2 


are : 


Res 


I  .  d 

<zp  +  d* 

l^p-0.0 

- (ip  -  e  T/2)  (z  -  Zp) 

~T  »  2e~T/2 
-2T 


(176) 


-0.0 


Res 


d 

(zp  +  l)z 

1  ,  -T  -  -T/2.  2  .  ,,  _2T  .  ,  -3T/2. 

*-(e  +  2e  )z  +  (3e  LL  +  4e  )z 

-T/2  "  dZP 

(z  )2(z  -  zz) 

1  -2T,  -T\  2 

e  (s  -  e  ) 

Ve 

L  P  p 

Ive- 

(177) 


Combining  Eqs.  176  and  177  gives  the  low-rate  discrete  transfer  function 
in  Eq.  178. 


104 


(1  -  e'1^)  Res 


Up-0.0 


-  e~T/2)2  z  +  (e~T  +  2e~T/2) 


(z  -  e-T)2 


,  *  -  e- 


O.OKz  +  2.61) 
(z  -  0.81) 2 


(178) 


Equation  174  is  now  solved  using  the  TXCONV  computer  program.  The 
required  input  data  are  outlined  below: 


Data  Card 


Data  Items 


.01*1,3 

.1053605155,-210721031 


-1.0, 0.0, 0.0, 0.0 

.9,0.0, .9,0.0 


The  output  for  this  example  is  shown  in  Fig.  26.  Both  the  high-rate  and 
low-rate  z-plane  transfer  functions  are  printed. 

A  second  example  will  consider  the  w'-plane  high-rate  to  low-rate 
transform  conversion  in  Ea.  179. 


(179) 


where 


[KG(w')T/3]- 


K  =  .01274481960 


T  =  .04  (sec) 
T/3  =  .12  (sec) 


The  high-rate  zeros  and  poles  for  the  w'-plane  transfer  function 
G(w')1^  are  given  by: 
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Zeros  (real,imag) 


Poles  (real,imag) 


(-.9999481996  , 
(.001830897924, 
(-.5657507592  , 
(-.5657507592  , 
(-4.956749255  , 
(-9.892804162  , 
(-14.00440188  , 
(-15.45466095  , 
(50.000000000  , 


0.0  ) 

0.0  ) 

6.861565379  ) 
-6.861565379) 


0.0  ) 

o.o  ) 

0.0  ) 

0.0  ) 

0.0  ) 


(-.5087865665  , 
(-.5087865665  , 
(-.001726986844, 
(-2.169528987  , 
(-2.169528987  , 
(-12.82548298  , 
(-16.86376037  , 
(-16.86376037  , 
(-10.52412126  , 


.3042290447  ) 
-.3042290447) 
0.0  ) 
3.364297796  ) 
-3.364297796) 
0.0  ) 
-12.01331830) 
12.01331830  ) 
0.0  ) 


A  listing  of  the  input  data  for  this  example  is  outlined  below: 


Data  Cards 

1  .01274481960,9,9 

2  .04,. 12 

3  WP 


4  Zeros  (real,imag) 

nth  Poles  (real,imag) 


The  output  for  this  w'-plane  example  is  shown  in  Fig.  27.  Figure  27a 
contains  the  high-rate  w'-plane  transfer  function  input  and  Fig.  27b  the 
low-rate  w'-plane  transfer  function  output. 

The  high-rate  w'  transfer  function  KG(w')^^  was  obtained  from  the 
DISCRET  computer  program  (Section  IV)  using  the  sampled  continuous 
system  defined  in  Eq.  180. 

KG(v')T/3  =  [k1G1(s)M3(s)]T/3  (1.80) 
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In  Eq.  180, 


M3  (s  ) 


1  _  e-sT/3 
s 


KX  =  -.6483736462  , 


and  the  s-plane  zeros  and  poles  for  G^(s)  are  given  by: 


_ Zeros  (real,imag) _ 

(-1.000000000  ,  0.0  ) 

(.001830897332,  0.0  ) 

(-5 .000000000  ,0.0  ) 

(-.5008733927  ,  6.832938756  ) 
(-.5008733927  ,  -6.832938756) 
(-15.00000000  ,  0.0  ) 

(-15.00000000  ,  0.0  ) 

(-10.000000000  ,  0.0  ) 


_ Poles  (real,imag) _ 

(-.5087852889  ,  .3042567928  ) 
(-.5087852889  ,  -.3042567928) 
(-.001726986844,0.0  ) 
(-2. Ibl077353  ,  3.365523190  ) 
(-2.161077353  ,  -3.365523190) 
(-13.11843178  ,  0.0  ) 
(-16.40426112  ,  -13. l3y42/24) 
(-16.40426112  ,  13.13942724  ) 
(-10.68380409  ,  0.0  ) 


These  parameters  were  used  with  the  WP  option  in  D1SCRET  to  obtain  the 
w'-plane  transfer  function  KG(w')T/^. 
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Figure  26.  T7C0NV  z -Plane  Output  for  Example 
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Figure  27s..  TXCONV  w'-Plane  Output  for  Example 


SECTION  VI 


SUMMARY 


The  high-speed  timeshared  digital  computer  with  a  large  memory  has 
made  it  possible  for  the  analyst  and  the  computing  machine  to  be  coupled 
closely  together  to  solve  a  large  variety  of  engineering  problems.  In 
system  analysis  and  design,  one  of  the  essential  ingredients  for  this 
coupling  is  a  library  of  properly  structured  analysis  and  synthesis  pro¬ 
grams  residing  in  the  machine.  Moreover,  easily  applied  and  accurate 
computer  programs  that  deal  with  discrete  or  hybrid  systems  are  becoming 
more  essential  as  a  result  of  the  increasing  use  of  digital  controllers 
In  automatic  control  systems.  Practical  hybrid  systems  contain  continu¬ 
ous  or  analog  elements  (e.g.,  plant,  process,  or  controlled  element) 
which  can  be  described  by  or  approximated  with  linear  differential  equa¬ 
tions  and  discrete  elements  (e.g.,  digital  controller)  which  are  inher¬ 
ently  defined  by  difference  equations.  Thus,  not  only  is  the  high¬ 
speed,  large  memory  digital  computer  used  to  analyze  and  design  discrete 
or  hybrid  systems,  but  the  small-scale  digital  computer  is  becoming  an 
integral  part  of  the  control  system  itself. 

The  fundamental  first  step  in  the  analysis  or  design  of  a  hybrid 
control  system  is  to  acquire  or  formulate  a  valid  linear  model  of  the 
system  being  considered.  This  includes  the  discretization  of  the  con¬ 
tinuous  elements  In  the  system  into  a  valid  discrete  domain.  The 
DISCRET  computer  program  presented  in  this  report  provides  this  discre¬ 
tization.  DISCRET  takes  a  continuous  element  expressed  as  an  s-plane 
transfer  function  and  transforms  it  into  the  z-,  w,  or  w'-plane  (see 
Eq.  58).  The  resulting  discrete  transfer  function  is  exact  as  opposed 
to  the  approximate  discrete  transfer  function  obtained  from  a  substitu— 
tion-for-s  approach  such  as  the  Tustin  or  first-difference  transforms. 
The  discrete  transfer  function  generated  by  DISCRET  defines  the  continu¬ 
ous  variables  (associated  with  the  continuous  element)  at  each  sampling 
Instant  of  the  sampler  device  in  the  system.  It  is  assumed  that  the 
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sampler  passes  the  continuous  signal  at  discrete  instances  equally 
spaced  in  time.  The  time  interval  T  between  these  samples  is  called  the 
sampling  period.  Once  calculated,  this  discrete  model  of  the  continuous 
element  can  be  readily  combined  with  the  inherent  discrete  elements  in 
the  system  (such  as  a  digital  controller)  to  provide  a  unified,  concise 
description  of  the  total  system.  An  inherently  discrete  element  in  a 
hybrid  system  is  first  modeled  with  a  recursion  or  difference  equation 
and  then  directly  converted  to  the  z-,  w-,  or  w'-plane  by  substituting 
the  z-n  delay  operator  for  each  discrete  term  in  the  difference  equa¬ 
tion.  Therefore,  all  the  elements  in  a  hybrid  system  can  be  described 
by  transfer  functions  in  the  z- ,  w~,  or  w'-plane.  Consequently,  similar 
frequency  and  time  domain  techniques  used  for  continuous  s-plane  systems 
can  be  applied  by  the  control  system  engineer  to  a  wide  variety  of  prac¬ 
tical  hybrid  systems. 

Many  practical  hybrid  control  systems  are  also  multi-rate  in  nature. 
That  is,  they  contain  samplers  which  operate  at  different  sample  rates. 
This  adds  additional  complexity  to  Che  analysis  and  design  procedures, 
but  concise,  definitive  methods  are  available  for  handling  multi-rate 
systems  (e.g..  Refs.  1-3).  However,  a  troublesome  situation  arises  when 
the  output  of  a  multi-rate  system  is  sampled  at  a  lower  rate  than  its 
input.  This  sampling  configuration  invariably  leads  to  the  requirement 
of  computing  a  low-rate  discrete  transform  from  a  given  high-rate  dis-. 
c.rete  transform  (see  Eqs.  76-78).  The  TXCONV  computer  program  presented 
in  this  report  calculates  this  complex  transformation  and  makes  it  a 
routine  operation  in  the  analysis  and  design  process.  Specif icaliy, 
TXCONV  takes  a  high-rate  discrete  transfer  function  expressed  in  the  z-, 
w-,  or  w'-plane  and  transforms  it  into  a  desired  low-rate  discrete 
transfer  function  in  the  z-,  w-,  or  w'-plane.  This  relieves  the  analyst 
from  the  computationally  involved  procedure  of  calculating  the  residues 
of  a  complex  integral. 

Both  the  DISCRET  and  -TXCONV  computer  programs  represent  essential 
tools  for  dealing  with  multi-rate,  hybrid  control  systems.  They  provide 
two  of  the  primary  techniques  used  to  formulate  a  unified  discrete  model 
of  a  complex,  hybrid  system.  However,  they  are  by  no  means  the  only 
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tools  or  techniques  required  by  the  analyst.  Other  computational  opera¬ 
tions  that  are  required  include  polynomial  and  transfer  function  manipu¬ 
lation  routines.  These  routines  are  needed  to  handle  the  individual 
discrete  transfer  functions  present  in  single-loop  and  multiloop  hybrid 
systems.  Analysis  routines  that  calculate  root  locus,  frequency  re¬ 
sponse,  and  time  response  in  the  discrete  domain  for  single-rate  and 
multi-rate,  hybrid  systems  are  also  needed.  Most  of  these  routines  for 
single-rate  discrete  or  hybrid  systems  are  already  available  in  the 
TOTAL  computer  program  (Ref.  13).  Additional  software  routines  that 
specifically  address  multi-rate  systems  and  the  special  features  in¬ 
volved  in  calculating  the  continuous  frequency  response  of  a  discretely 
excited  system  (Ref.  1)  are  presently  under  development  at  AFWAL/FIGC. 

Both  DISCRET  and  TXCONV  have  been  integrated  into  the  total  computer 
program.  DISCRET  and  TXCONV  are  two  of  the  over  100  options  available 
during  the  interactive  execution  of  the  TOTAL  program.  This  interactive 
feature  allows  a  close  coupling  between  the  analyst  and  the  computing 
machine  (digital  computer)  such  that  a  real-time  dialog  between  the  two 
can  be  effectively  carried  out.  This  results  in  a  more  effective  and 
efficient  visage  of  the  computing  machine  and  improves  the  accuracy  and 
speed  of  the  analysis  and  synthesis  process. 
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